2017-07-19 25 views
0

我创建了一些MatLab代码,它使用两个给出相同平面波的不同表达式绘制平面波。第一个表达式是在笛卡尔坐标中,并且正常工作。然而,第二个表达式是在极坐标中,当我在这种情况下计算平面波时,图是失真的。这两块地块看起来应该是一样的。那么,我在做极地坐标转换时做错了什么?MatLab中有两个相同的波形图,但转换为极坐标后创建的图形会变形吗?

function Plot_Plane_wave() 
clc 
clear all 
close all 

%% Step 0. Input paramaters and derived parameters. 
alpha = 0*pi/4;    % angle of incidence 
k = 1;      % wavenumber 
wavelength = 2*pi/k;  % wavelength 


%% Step 1. Define various equivalent versions of the incident wave. 
f_u_inc_1 = @(alpha,x,y) exp(1i*k*(x*cos(alpha)+y*sin(alpha))); 
f_u_inc_2 = @(alpha,r,theta) exp(1i*k*r*cos(theta-alpha)); 


%% Step 2. Evaluate the incident wave on a grid. 
% Grid for field 
gridMax = 10; 
gridN = 2^3; 
g1 = linspace(-gridMax, gridMax, gridN); 
g2 = g1; 
[x,y] = meshgrid(g1, g2); 

[theta,r] = cart2pol(x,y); 

u_inc_1 = f_u_inc_1(alpha,x,y); 

u_inc_2 = 0*x; 
for ir=1:gridN 
    rVal = r(ir); 
    for itheta=1:gridN 
     thetaVal = theta(itheta); 
     u_inc_2(ir,itheta) = f_u_inc_2(alpha,rVal,thetaVal); 
    end 
end 

%% Step 3. Plot the incident wave. 
figure(1); 
subplot(2,2,1) 
imagesc(g1(1,:), g1(1,:), real(u_inc_1)); 
hGCA = gca; set(hGCA, 'YDir', 'normal'); 
subplot(2,2,2) 
imagesc(g1(1,:), g1(1,:), real(u_inc_2)); 
hGCA = gca; set(hGCA, 'YDir', 'normal'); 

end 

Plots expected to be the same

回答

1

你的错误就在于,你的循环只有通过rtheta第一gridN值去。相反,您需要逐步筛选ixiy的索引,并提取矩阵的rValthetaVal以及rtheta

您可以将循环改为

for ix=1:gridN 
    for iy=1:gridN 
     rVal = r(ix,iy);   % Was equivalent to r(ix) outside inner loop 
     thetaVal = theta(ix,iy); % Was equivalent to theta(iy) 
     u_inc_2(ix,iy) = f_u_inc_2(alpha,rVal,thetaVal); 
    end 
end 

这给预期的图表。

Plots from corrected code


或者,你可以在你的内联函数喂养矩阵简化代码。要做到这一点,你将不得不在f_u_inc_2使用的元素单元的产品,而不是.*矩阵乘法*的:

alpha = 0*pi/4; 
k = 1; 
wavelength = 2*pi/k; 

f_1 = @(alpha,x,y) exp(1i*k*(x*cos(alpha)+y*sin(alpha))); 
f_2 = @(alpha,r,theta) exp(1i*k*r.*cos(theta-alpha)); 
% Change       v 
f_old = @(alpha,r,theta) exp(1i*k*r *cos(theta-alpha)); 

gridMax = 10; 
gridN = 2^3; 

[x,y] = meshgrid(linspace(-gridMax, gridMax, gridN)); 
[theta,r] = cart2pol(x,y); 

subplot(1,3,1) 
contourf(x,y,real(f_1(alpha,x,y))); 
title 'Cartesian' 
subplot(1,3,2) 
contourf(x,y,real(f_2(alpha,r,theta))); 
title 'Polar' 
subplot(1,3,3) 
contourf(x,y,real(f_old(alpha,r,theta))); 
title 'Wrong' 

Output image from code