吸声材料仿真实验中的角度问题

李嘎嘎

clear all
f=1000;
w=2*pi*f;
c=344;
k0=w/c;
rho=1.29;%空气密度
Z0=rho*c;%空气特征阻抗
delta=100;%单位：cgs（流动阻抗）
%%
x1=0;
y1=0;
%%
dxy=0.05;%网格间距；
N=5.5;%网格数目；
Lx=N*dxy;
Ly=N*dxy;
x=-Lx:dxy:Lx;
y=-Ly:dxy:Ly;
[X,Y]=meshgrid(x,y);
zr=0;%材料表面（默认为Z=0的表面）
zh=0.0;%测量面
zs=0.05;%声源面
r=sqrt((X-x1).^2+(Y-y1).^2);%xy方向的位置
R1=sqrt(r.^2+(zs-zh).^2);%zr、zh为重建与测量面在Z方向的位置
%吸声材料基本参数设定
Z1=Z0*(1+9.08*(f/delta)^(-0.75)+1i*11.9*(f/delta)^(-0.73));
k1=k0*(1+10.8*(f/delta)^(-0.70)+1i*10.3*(f/delta)^(-0.59));

%% 就是这两个角度,让我发狂!!
THETA=asin(sqrt(X.^2+Y.^2)./R1);
PHI=atan(Y./X);
%% 就是这两个角度,让我发狂!!


kx=k0.*sin(THETA).*cos(PHI);
ky=k0.*sin(THETA).*sin(PHI);
% kr=sqrt(kx.^2+ky.^2);
% kr=k0.*sin(theta);
Cr_mol=Z1*cos(THETA)-Z0*sqrt(1-sqrt(k0/k1)*(sin(THETA)).^2);
Cr_dem=Z1*cos(THETA)+Z0*sqrt(1-sqrt(k0/k1)*(sin(THETA)).^2);
Cr_kr=Cr_mol./Cr_dem;
P_in=exp(1i*k0*R1)./(4*pi*R1);
%
kr=eps:0.5:100;
r=reshape(r,[],1);
Cr_kr=reshape(Cr_kr,[],1);
for jj=1:length(r)
    kz=sqrt(k0^2-kr.^2);
    for mm=1:length(kz)
        GREEN=exp(1i*kz(mm)*(zs+zh))/(4*pi*kz(mm));
        F_x=1i*GREEN*Cr_kr(jj).*bessely(0,r(jj)*kr).*kr;
        P_re(jj)=trapz(kr,F_x);
    end
end
%
P_RE=reshape(P_re,length(x),length(y));
PH=P_in+P_RE;
figure(1)
subplot(1,2,1),surf(abs(P_RE))
title('反射')
subplot(1,2,2),surf(abs(P_in))
title('入射')
现在的问题是程序是没有问题,但是中间的逻辑我有点搞不清了,正如我前一个帖子所述:到底这几个角度是怎么回事?
理论来自于:Spatial Fourier-transform method for measuring re...
问题求助:关于参数THETA的理论意义,在论文里的解释就是K空间里波数K0相对于x,y,z方向上的一个夹角,可是在仿真实验结果里却存在一个,声波相对于材料表面的入射角问题,我就找不到这个所谓的入射角到底是他妈的哪个角啦!
艹!!!急死我了!!!