matlab shading 网格密度,matlab关于shading interp的问题

%**********************************************************************

% Declaring Contants

%**********************************************************************

tic

c=3e8; % Speed of light 光速

f=1e9; % Operating frequency 频率

Lamda=c/f; % Wavelength

ko=(2*pi)/Lamda; % Wavenumber 波数

S=5000; % Ground conductivity in mS/m 磁导率

Er=80; % Relative dielectric constant 相对介电常数

Dz=Lamda; % Delta Z in meter z方向间隔等于波长

Dy=Lamda; % Delta Y in meter y方向波长

sigmaz=Lamda; % Current source standard deviation 标准差等于波长

%**********************************************************************

% Input Parameters

%**********************************************************************

Ht=10; % Transmitting antenna height 天线高度

D=1000; % Range from the transmitter x方向传播距离

Dx=100; % Incremental range (Dx) in meter 步长

Ny=512; % Sample size in the y-direction 虚列

Nz=512; % Sample size in the z-direction

%**********************************************************************

% Define range in y-direction and z-direction

%**********************************************************************

ymax=Ny*Dy/2; % maximum range in y-direction y轴最大值

y1=0:Dy:ymax; % First half of y range incrementally increase y正向增量

y=[y1 zeros(1,Ny/2-1)]; % y array

zmax = Nz*Dz/2; % Find maximum z range (Zmax)

z1=0:Dz:zmax; % First half of the range of Zmax+1

z=[z1 zeros(1,Nz/2-1)]; % Construct a full z array

%**********************************************************************

% Define Wavenumbers in Spatial Domain

%**********************************************************************

kymax=pi/Dy; % Maximum wavenumber in y-direction 频域中ky最大值

Dky=2*kymax/Ny; % detla ky  ky间隔

ky1=-kymax:Dky:kymax; % Range of ky

ky=[ky1(:,(Ny/2)+1:Ny+1) ky1(:,2:Ny/2)];

kzmax=pi/Dz; % Maximum wavenumber in z-direction

Dkz=2*kzmax/Nz; % detla kz

kz1=-kzmax:Dkz:kzmax; % Range of kz

kz=[kz1(:,(Nz/2)+1:Nz+1) kz1(:,2:Nz/2)];

%**********************************************************************

% Compute Wavenumbers in Frequency Domain, kx

%**********************************************************************

Ky=meshgrid(ky,1:Nz); % Create a NyxNz matrix of ky row-repeat

Kz=meshgrid(kz,1:Ny).'; % Create a NyxNz matrix of kz columnrepeat

kx=sqrt(ko^2-Ky.^2-Kz.^2); % Compute theNy x Nz kx matrix

clear Ky; % Clear matrices Ky and Kz from memory

clear Kz; % To make the algorithm run faster

%**********************************************************************

% Compute reflection coefficient

%**********************************************************************

Erc=Er+i*18*S/(f/1e6); % Complex dielectric 介电常数

Zs=1/sqrt(Erc); % Impedance 阻抗

Gamma=(kz-ko*Zs)./(kz+ko*Zs); % Complex reflection coefficient 反射系数

%gamma=-ones(Nz,Ny); % Setting gamma = -1

gamma=meshgrid(Gamma,1:Ny).'; % Create reflection coefficient matrix

% by column repeating

%**********************************************************************

% Define the Hanning Window

%**********************************************************************

Hya=[]; % An empty vector

for t=0:Ny/2; % Define Ny/2 + 1 points

if (t>=0 & t<=3*Ny/8) % Define first 3Ny/8 +1 points

hy=1;

elseif (t>=3*Ny/8 & t<=Ny/2)% Define next Ny/8 points

hy=(sin(4*pi*t/Ny))^2;

end

Hya=[Hya hy]; % Construct the first half of Hanning

end % window of 1+Ny/2 elements

Yy=fliplr(Hya(:,2:Ny/2)); % Flip left to right of the Hanning

% Ny/2 element for total of [1x(Ny/2-1)]

HY=[Hya Yy]; % Hanning window in y-direction

Hzb=[]; % An empty vector

for t=0:Nz/2; % Define Nz/2 + 1 points

if (t>=0 & t<=3*Nz/8) % Define first 3Nz/8 +1 points

hz=1;

elseif (t>=3*Nz/8 & t<=Nz/2)% Define next Nz/8 points

hz=(sin(4*pi*t/Nz))^2;

end

Hzb=[Hzb hz]; % Construct the first half of Hanning

end % window of 1+Nz/2 elements

Yz=fliplr(Hzb(:,2:Nz/2)); % Flip left to right of the Hanning

% windows from the second element to

% Nz/2 element for total of [1x(Nz/2-1)]

HZ=[Hzb Yz]'; % Hanning window in z-direction

Hmy=meshgrid(HY,1:Nz); % Row repeat

Hmz=meshgrid(HZ,1:Ny)'; % Column repeat

H3D=(Hmy.*Hmz); % 3D Hanning window

%**********************************************************************

% Initial H-field at x=0 初始场

%**********************************************************************

g_tilda=exp(-kz.^2*sigmaz^2/2); % initial g tilda

hye0_tilda=g_tilda.*cos(kz*Ht); % Initial hye_tilda(0+,ky,kz) field

hyo0_tilda=-i*g_tilda.*sin(kz*Ht); % Initial hyo_tilda(0+,ky,kz) field

Hye0_tilda=meshgrid(hye0_tilda,1:Ny).';% Initial even H-field, column repeat

Hyo0_tilda=meshgrid(hyo0_tilda,1:Ny).';% Initial odd H-field, column repeat

Hy0_tilda=0.5*(1+gamma).*Hye0_tilda+0.5*(1-gamma).*Hyo0_tilda;% Include the reflection coefficient

Hy_tilda=Hy0_tilda.*H3D;% Hy_tilda(0+,ky,kz) at x=0

Exp2=exp(i*kx*Dx); % Marching range 每循环一次都乘以它

%**********************************************************************

% 3D Parabolic Equation Basic Algorithm 方程循环计算

%**********************************************************************

x=0:Dx:D-Dx;

nn=length(x);

for k=1:nn

Hy_tilda_Dx=Hy_tilda.*Exp2; % Hy_tilda(Dx,ky,kz)

Hy=Dkz*Dky*(Ny*Nz*ifft2(Hy_tilda_Dx))/(2*pi)^2;% ifft2 wrt to ky and kz

Hy_H=Hy.*H3D; % Apply the Hamming window

Hh(k,:,:)=abs(Hy_H(0:512,0:512));

Hyz1=Hy_H(1:Nz/2+1,:); % H-field for z > 0

Hyz0=flipud(Hy_H(2:Nz/2,:)); % H-field for z < 0

Hye=[Hyz1; Hyz0]; % Even H-field

Hyo=[Hyz1; -Hyz0]; % Odd H-field

Hye_tilda=(Dz*Dy)*fft2(Hye); % Take the FFT of the even H~

Hye_tilda_H=Hye_tilda.*H3D; % Apply the Hanning window in freq

Hyo_tilda=(Dz*Dy)*fft2(Hyo); % Take the FFT of the odd H~

Hyo_tilda_H=Hyo_tilda.*H3D; % Apply the Hanning window in freq

% Multiply by the reflection coefficient

Hye_tilda_g=0.5*(1+gamma).*Hye_tilda_H;

% Multiply by the reflection coefficient

Hyo_tilda_g=0.5*(1-gamma).*Hyo_tilda_H;

Hy_tilda1=Hye_tilda_g + Hyo_tilda_g;% Hy_tilda(x,ky,kz)

Hy_tilda=Hy_tilda1.*H3D; % Apply the Hanning window

Hh(k,:,:)=abs(ifft(Hy_tilda));

end

slice(Hh,100,5,5);

shading interp;colormap hsv

我这程序画出的三维图总感觉，xy轴反了！本是想画电波源在原点位置的发射效果图，但是图出来就感觉xy轴换了，请求高手解答。