INTEGARTION OF BESSELH function

Hi

i make a code where i integrate the hankel function besselh(0, 2)

function z = Escat(r,phi)

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

[Is]=currentMoM()

Phim=zeros(N+1);

FINAL=0;

for jj=1:N

%Phi0(jj)=(jj-1).*(pi./N);

Phi0(jj)=(jj-1).*(2.*pi./N);

FINAL=FINAL-(k0.*Z0./4).*Is(jj).*ra.*integral(@(xx)besselh(0,2,k0.*sqrt(ra.^2+r.^2-2.*r.*ra.*(phi-xx))),Phi0(jj),2*pi/N+Phi0(jj));

end

z=FINAL;

end

when i plot the function in the range φ=0 : 2π with r=const the valus of Escat are extremely bih value ...

clear all

clc

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

%ph_i=pi/2;

rho=10.*ra ;

phi=0:pi/180:2*pi;

Es=zeros(length(phi));

%phi=0:pi/200:pi/3

for jj=1:length(phi)

Es(jj)=abs(Escat(rho,phi(jj)));

end

plot(phi*(180./pi),Es,'b--')

hold on

%plot(phi*(180./pi),abs(integ),'b--')

plot(phi,abs(Escattheory_new(rho,phi)),'r-')

xlabel('$\phi$','Interpreter','latex')

ylabel('$|E_{s}|$','Interpreter','latex' )

%legend('MoM', 'Theory')

hold off

can i check if the the integration is ok ?

the parameter fucntion is

function [f,N,ra,k0,Z0] = parameter()

%UNTITLED Summary of this function goes here

c0=3e8;

Z0=120.*pi;

ra=1;

N=80;

f=300e6;

lambda=c0./f;

k0=2*pi./lambda;

end

thank you

[Steven Lord]

You haven't provided two of the functions used by your code: currentMoM and Escattheory_new.

You wrote "the valus of Escat are extremely bih value ..." -- I assume that's a typo and you meant to say "big value". How big is "extremely big" and how big were you expecting the values to be? Why did you expect the values to be so different from the values you actually received?

[george veropoulos]

thank for your reply

i want to compare the function Escat with Escattheory_new

function y=Escattheory_new(r,phi)

%UNTITLED Summary of this function goes here

% Detailed explanation goes here

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

E0=1;

Z1=k0.*r;

Z2=k0.*ra;

final=0 ;

for kk=0:200

final=final-E0*(besselj(kk,Z2,1).*besselh(kk,2,Z1)./besselh(kk,2,Z2))*(-1j).^kk.*cos(kk.*phi)*e_n(kk);

%z=vpa(symsum(factor,k,0,inf),3)

end

y=final;

end

end the values ar not the same also Near the φ=2*π the Escat take high values 10^16

George

[george veropoulos]

and the cuurentMoM

function [Is]=currentMoM()

%UNTITLED2 Summary of this function goes here

% Detailed explanation goes here

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

gamma_const=1.781;

Phi0=zeros(N);

e=exp(1);

dfpm=2.*pi./N;

for jj = 1:N

Phi0(jj)=(jj-1).*(2.*pi./N);

end

% delta_c(i) = sqrt((pos(i,1) - pos(i+1,1))^2 + (pos(i,2) - pos(i+1,2))^2);

lmn = zeros(N);

%zmn = zeros(N);

gm = zeros(1,N);

zmn = zeros(N);

%vim = zeros(1,N);

%vsn = zeros(1,N);

coeif=(Z0.*k0./4).*ra.*dfpm;

coeifn=(Z0./2).*sin(k0.*ra.*dfpm./2);

for index_i = 1:N

for index_j = 1:N

if index_i == index_j

lmn(index_i,index_j) =coeif.*(1 - 1i.*(2./pi).*(log((gamma_const.*k0.*ra.*dfpm)./4)-1));

else

R(index_i,index_j)=ra.*sqrt(2.-2.*cos(Phi0(index_j)-Phi0(index_i)));

lmn(index_i,index_j) =coeifn.*besselh(0,2,k0.*R(index_i,index_j));

end

zmn(index_i,index_j) = lmn(index_i,index_j);

gm(index_i) =Efieldin(Phi0(index_i));

end

%vim(index_i) = delta_c(index_i) * exp(j*k0*(xm(index_i)*cos(phi_i)+ym(index_i)*sin(phi_i)));

%vsn(index_i) = delta_c(index_i) * exp(j*k0*(xm(index_i)*cos(phi_s)+ym(index_i)*sin(phi_s)));

end

W = linsolve(zmn,gm');

for ii=1:N

Is(ii)=W(ii);

end

y= Is

end

[george veropoulos]

function y=e_n(k)

if k==0

y=1;

elseif k~=0

y=2;

end

end

[Steven Lord]

If you're adding information, please put it as a comment rather than a separate answer. Leave the answer for posts that may be a solution to the problem, to make them easier to distinguish from commentary.