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

[Walter Roberson]

Note that function e_n leaves y undefined in the case where k is NaN.

If k cannot be NaN then there is no point in using elseif -- just use else

Answer 1

[David Goodmanson]

Hi george,

I don't know what the values of your parameters are, but it appears that the problem with Escat is that

besselh(0,2,k0.*sqrt(ra.^2+r.^2-2.*r.*ra.*(phi-xx)))  

is missing a cosine factor and should be

besselh(0,2,k0.*sqrt(ra.^2+r.^2-2.*r.*ra.*cos(phi-xx)))

Without the cosine, it's possible for the argument of the sqrt to become negative, which causes the argument of besselh to become imaginary, which leads to very large values in the result.

[george veropoulos]

thank you very much!!! the cos was the problem