h=0.0007;
c=3e8;%光速
er=2.33;%相对介电常数
e0=8.8542*10^(-12);%真空介电常数
f=2e9;%frequency
k0=2*pi*f/c;%自由空间波数
probex=x1;probey=y1;probez=z1;%监测点坐标
dipolex=x2;dipoley=y2;dipolez=z2;%等效偶极子坐标
for m=1:882;
for n=1:128;
%求柱坐标系
fi(m,n)=atan((probey(m)-dipoley(n))/(probex(m)-dipolex(n)));
r(m,n)=sqrt((probex(m)-dipolex(n))^2+(probey(m)-dipoley(n))^2);
z(m,n)=probez(m)-dipolez(n);
%转移系数矩阵
if mod(m,2)==1
if mod(n,2)==1
fun1=@(x)(x.^3./((sqrt(x.^2-k0^2)+sqrt(x.^2-er*k0^2).*coth(h*sqrt(x.^2-er*k0^2))).*(er*sqrt(x.^2-k0^2)+sqrt(x.^2-er*k0^2).*tanh(h*sqrt(x.^2-er*k0^2))))).*exp(-z(m,n)*sqrt(x.^2-k0^2)).*(besselj(0,x*r(m,n))-(besselj(1,x*r(m,n))./(x*r(m,n))));
integral11(m,n)=quadgk(fun1,0,1.0001*k0);
integral21(m,n)=integral11(m,n);
integral12(m,n)=quadgk(fun1,1.000145*k0,inf);
integral22(m,n)=integral12(m,n);
a(m,n)=(((er-1)*sin(2*fi(m,n)))/(2*pi))*(integral11(m,n)+integral12(m,n));
else
fun21=@(x)besselj(0,x*r(m,n)).*x.*sqrt(x.^2-k0^2).*exp(-z(m,n)*sqrt(x.^2-k0^2))./(sqrt(x.^2-k0^2)+sqrt(x.^2-er*k0^2).*coth(h*sqrt(x.^2-er*k0^2)));
integral11(m,n)=quadgk(fun21,0,1.0001*k0);
integral12(m,n)=quadgk(fun21,1.000145*k0,inf);
fun22=@(x)(x.^3./((sqrt(x.^2-k0^2)+sqrt(x.^2-er*k0^2).*coth(h*sqrt(x.^2-er*k0^2))).*(er*sqrt(x.^2-k0^2)+sqrt(x.^2-er*k0^2).*tanh(h*sqrt(x.^2-er*k0^2))))).*exp(-z(m,n)*sqrt(x.^2-k0^2)).*((sin(fi(m,n)))^2*besselj(0,x*r(m,n))+cos(2*fi(m,n))*(besselj(1,x*r(m,n))./(x*r(m,n))));
integral21(m,n)=quadgk(fun22,0,1.0001*k0);
integral22(m,n)=quadgk(fun22,1.000145*k0,inf);
a(m,n)=(integral11(m,n)+integral12(m,n))/(2*pi)-(((er-1)*sin(2*fi(m,n)))/(2*pi))*(integral21(m,n)+integral22(m,n));
end
else
if mod(n,2)==1
fun31=@(x)besselj(0,x*r(m,n)).*x.*sqrt(x.^2-k0^2).*exp(-z(m,n)*sqrt(x.^2-k0^2))./(sqrt(x.^2-k0^2)+sqrt(x.^2-er*k0^2).*coth(h*sqrt(x.^2-er*k0^2)));
integral11(m,n)=quadgk(fun31,0,1.0001*k0);
integral12(m,n)=quadgk(fun31,1.000145*k0,inf);
fun32=@(x)(x.^3./((sqrt(x.^2-k0^2)+sqrt(x.^2-er*k0^2).*coth(h*sqrt(x.^2-er*k0^2))).*(er*sqrt(x.^2-k0^2)+sqrt(x.^2-er*k0^2).*tanh(h*sqrt(x.^2-er*k0^2))))).*exp(-z(m,n)*sqrt(x.^2-k0^2)).*((cos(fi(m,n)))^2*besselj(0,x*r(m,n))-cos(2*fi(m,n))*(besselj(1,x*r(m,n))./(x*r(m,n))));
integral21(m,n)=quadgk(fun32,0,1.0001*k0);
integral22(m,n)=quadgk(fun32,1.000145*k0,inf);
a(m,n)=(integral11(m,n)+integral12(m,n))/(2*pi)-(((er-1)*sin(2*fi(m,n)))/(2*pi))*(integral21(m,n)+integral22(m,n));
else
fun4=@(x)(x.^3./((sqrt(x.^2-k0^2)+sqrt(x.^2-er*k0^2).*coth(h*sqrt(x.^2-er*k0^2))).*(er*sqrt(x.^2-k0^2)+sqrt(x.^2-er*k0^2).*tanh(h*sqrt(x.^2-er*k0^2))))).*exp(-z(m,n)*sqrt(x.^2-k0^2)).*(besselj(0,x*r(m,n))-(besselj(1,x*r(m,n))./(x*r(m,n))));
integral11(m,n)=quadgk(fun4,0,1.0001*k0);
integral21(m,n)=integral11(m,n);
integral12(m,n)=quadgk(fun4,1.000145*k0,inf);
integral22(m,n)=integral12(m,n);
a(m,n)=(((er-1)*sin(2*fi(m,n)))/(2*pi))*(integral11(m,n)+integral12(m,n));
end
end
end
end