MATLAB преминава от Cart към Pol обратно към Cart coords за цилиндричен график

1. Какво ИСКАМ да направя:
(i) Използвайте вход n, за да генерирате n*n декартова мрежа

[x y] = meshgrid(linspace(-1,1,n));  

(ii) Генериране на полярни координати

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

(iii) Оценява функция в цилиндрични координати

z = f(theta,r);  

(iv) Начертайте резултата, като използвате (да кажем) pcolor (или сърфиране, или нещо друго)

pcolor(x,y,abs(z).^2) %The function is complex, a Laguerre-Gauss to be exact.  

2. Какво МОГА да направя... Единственият начин да накарам диаграмите да работят е като започна с моите полярни параметри и оттам се върна към декартово:
(i)
Дефиниране на параметри

r=linspace(0,1,n); theta=linspace(0,2*pi,n);  

(ii) Създайте и двете мрежи и оценете f

[theta r]=meshgrid(theta,r);  
[x y]=pol2cart(theta,r);  
z=f(theta,r);  

(iii) Парцел

pcolor(x,y,abs(z).^2)  

ПРОБЛЕМЪТ е, че сега моята решетка е кръгла и бих искал да оценя функцията навсякъде ВЪРХУ ПРАВОЪГЪЛНА решетка (защото моят анализ зависи от наличието на квадратни пикселни масиви). Повтарянето, използвайки метод 2 по-горе, получавам кръгова графика, описана в квадрат; представете си черен кръг с бяло по краищата... но ИСКАМ да оценя функцията в тази "бяла" област. ОБАЧЕ използването на метод 1 НЕ работи -- цялата функция е объркана, когато чертая (Просто потърсете в Google Laguerre-Gauss modes, за да видите как трябва да изглеждат графиките).

Искам да мога да започна с правоъгълна мрежа и да присвоя на всяка точка полярна координата, вместо да започна с полярни координати и да им присвоя всички декартови точки.

Бъркам се с това на разстояние от дълго време и не мога да разбера как да заобиколя този на пръв поглед прост проблем.

Редактиране 1
Изглежда, че проблемът е в начина, по който се генерират координатните матрици. По-долу съм публикувал екранни снимки на прост пример 3 на 3, илюстриращ как подход 1 и подход 2 генерират различни числа.

Как да направим тези номера съвместими?

Нямам точки за репутация, така че не мога да кача изображенията директно... връзките по-долу показват примера 3 на 3... вижте коментарите за връзки към действителни изображения на графиките на Лагер-Гаус, които се опитвам да направя...

приложи cart2pol
приложи pol2cart

Редактиране 2

В момента резултатът от подхода (1.) дава грешни резултати, както е показано тук:

желан подход, грешен резултат

Вторият подход дава правилните изображения, за съжаление това е само кръг, а не целият квадрат. Показано е тук:

приложен подход, ограничен резултат

3D диаграми и на двата подхода са показани тук - само цветната част на горната фигура е вярно.

Редактиране 3

Ето екранна снимка на функцията f, която се използва по-горе. Обърнете внимание, че изисква повече входни параметри, отколкото само r,theta. Типичните стойности са:

w0 = 0.5;
p = 0;
l = 5;

Функцията C дава нормализация, а L са полиноми на Лагер. И двете функции са щателно тествани и дават очакваните резултати.

Редактиране 4
Ето достатъчно код, за да стартирам изрично моя пример z=U(0,5,r,phi,w0)+U(0,-5,r,phi,w0);. Самият сюжет е даден от pcolor(x,y,abs(z).^2).

Имайте предвид, че функцията Lpl() е вмъкната като коментар. Това ще трябва да бъде запазено като собствен m-файл, за да може функцията U да работи правилно.

%% Laguerre-Gauss Modes U = U(p,l,r,phi,w0)
%  Source: OAM theory paper section 2.A eqn 1.
%  Assuming POLAR coordinates and evaluating AT beam waist.
% -- That is, z=0 for w(z)=w0(sqrt(1+z/zR)) 
% ---- ie, w(0) = w0
% Assuming z=0 also renders the Gouy phase arctan(z/zR) irrelevant.
% Note: Rayleigh Range zR is not explicitly defined because z=0 --> it is irrelevant too.
% Since zR is the only wavelength dependent term, wavelength also doesn't
% matter.

function out = U(p,l,r,phi,w0)
%Function handles for clarity
e = @(x) exp(x);
C = @(p,l) sqrt((2*factorial(p))/(pi*factorial(p+abs(l))));
L = @(p,l,z) Lpl(p,l,z);

%% Lpl() FUNCTION
% function out = Lpl(p,l,z)
% 
% l=abs(l);
% LL=0;
% for mm=1:p+1
%     m=mm-1;
%     L=LL;
%     LL= L+((-1)^m)*(factorial(p+l)/(factorial(p-m)*factorial(l+m)*factorial(m)))*(z.^m); 
% end
% out = LL;

%% 

out = (C(p,l)/w0)*...
    (((sqrt(2).*r)/w0)^abs(l))*...
    (e((-r.^2)/w0^2))*...
    (L(p,l,((2.*r.^2)/w0^2)))*...
    (e((-1)*1i*l.*phi)); ``

person caseyalan    schedule 21.05.2013    source източник
comment
Логично е, че числата няма да са напълно еднакви, но кръг с радиус 1 трябва лесно да се преобразува в декартова мрежа, която обхваща (-1,1) във всяка посока. В моите очи аз оценявам функцията в същия домейн (по същество) във всеки случай. Очевидно това не е вярно... Не виждам какво пропускам.   -  person caseyalan    schedule 22.05.2013
comment
Благодаря ви много за помощта!   -  person caseyalan    schedule 22.05.2013


Отговори (1)


Редактиране
Отговорът беше пренаписан въз основа на кода, предоставен в Редактиране 4 на въпроса.

Мисля, че проблемът произтича от функцията U. Не прилагате елементни операции към всички части на уравнението . Ако го промените на:

out = (C(p,l)./w0).* ...              % here it's a .* instead of *
    (((sqrt(2).*r)./w0).^abs(l)).* ... % here it's a .* instead of *
    (e((-r.^2)./w0.^2)).* ...          % here it's a .* instead of *
    (L(p,l,((2.*r.^2)./w0.^2))).* ...  % here it's a .* instead of *
    (e((-1)*1i*l.*phi));

Получавате следните два резултата, показани по-долу.

Тази фигура използва въведени декартови координати:

въведете описание на изображението тук

И тази фигура използва полярните координати:

въведете описание на изображението тук

„По-грубата“ разделителна способност на втората фигура се дължи на по-малко подходящата разделителна способност на решетката. Но по същество разрешавате същите функции.

person Schorsch    schedule 22.05.2013
comment
Това е!!! Тествах вашето решение върху около 30 комбинации и всички дават очаквани резултати. Това всъщност беше елементарна грешка при поелементно умножение. - person caseyalan; 22.05.2013
comment
ТХХХАННККК ВИУУУУУ!!! Моите летни изследвания най-накрая могат да започнат (тъй като тези функции са основата - буквално - за всичко, което правя) :) - person caseyalan; 22.05.2013
comment
Искам да съм сигурен, че правя всичко необходимо, за да ти осигуря максимален реквизит за това. Кажете ми, ако забравям нещо (ново в стека). - person caseyalan; 22.05.2013
comment
@caseyfitz: заповядайте! - след като приемете отговора, вие сте добре. - person Schorsch; 22.05.2013