2017-07-19 29 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