Edits by Melissa Guidry fixing Cartesian coordinatesHEADmaster

1 files changed, 33 insertions, 22 deletions

diff --git a/transverse/decompose.m b/transverse/decompose.m
index e9b583a..7b8fd53 100644
--- a/transverse/decompose.m
+++ b/transverse/decompose.m
@@ -1,9 +1,9 @@
 %----------------------------------------------------------------------------------------------

 % PROGRAM: decompose

-% AUTHOR:  Andri M. Gretarsson

+% AUTHOR:  Andri M. Gretarsson, edited by Melissa Guidry

 % DATE:    7/10/04

 %

-% SYNTAX: [coeffs,tmat]=decompose(z1,domain,type,terms,[q <,lambda,accuracy>]);

+% SYNTAX: [coeffs,tmat]=decompose(z1,domain,type,terms,[q <,lambda,accuracy>],coordtype);

 %           <...> indicates optional arguments

 %

 % Calculates the coefficients of the decomposition of the function z1 into

@@ -67,6 +67,11 @@
 %             of 1.54 would be rounded to 1.5 while a coefficient of 1.56 would be 

 %             rounded to 1.8.  Accuracy of 0 applies no rounding.  Specifying

 %             accuracy other than "0" does not speed up the calculation.

+% coordtype = (string) label for the type of coordinates supplied in r and theta. If 

+%             this argument is equal to 'pol', r and theta are assumed 

+%             to be polar coordinates (r,theta) as described above.  If on the other 

+%             hand coordtype='cart' then r is assumed to be the cartesian x coordinate 

+%             and theta the cartesian y coordinate.

 %

 % OUTPUT ARGUMENTS:

 % -----------------

@@ -88,7 +93,7 @@
 %       [x,y]=meshgrid([-pi/2:0.1:pi/2],[-pi/2:0.1:pi/2]); z1=cos(sqrt(x.^2+y.^2));

 %       clear domain; domain(:,:,1)=x; domain(:,:,2)=y;

 %       w=pi/4; R=1e3; lambda=1e-6; q=(1./R - i* lambda./pi./w.^2).^(-1);

-%       [coeffs,tmat]=decompose(z1,domain,'hg',140,[q,lambda,1e-6])

+%       [coeffs,tmat]=decompose(z1,domain,'hg',140,[q,lambda,1e-6],'cart')

 %       z1recomposed=recompose(domain,'hg',coeffs,[q,lambda,1e-6]);

 %       subplot(3,1,1); h=pcolor(x,y,abs(z1).^2); set(h,'EdgeColor','none'); axis square; colorbar

 %       subplot(3,1,2); h=pcolor(x,y,abs(z1recomposed).^2); set(h,'EdgeColor','none'); axis square; colorbar

@@ -96,9 +101,15 @@
 %

 % Last updated: July 18, 2004 by AMG

 %----------------------------------------------------------------------------------------------

-% SYNTAX: [coeffs,tmat]=decompose(z1,domain,type,terms,[q <,lambda,accuracy>]);

+% SYNTAX: [coeffs,tmat]=decompose(z1,domain,type,terms,[q <,lambda,accuracy>],coordtype);

 %----------------------------------------------------------------------------------------------

-function [coeffs,tmat]=decompose(z1,domain,type,terms,params)

+function [coeffs,tmat]=decompose(z1,domain,type,terms,params,coordtype)

+

+if strcmpi(coordtype,'cart')

+    metric = 1;

+elseif strcmpi(coordtype, 'pol')

+    metric = domain(:,:,1);

+end

 

 % HERMITE GAUSSIAN EXPANSION --------------------------------------------------------------------------

 if strcmpi(type,'hg')

@@ -109,10 +120,10 @@ if strcmpi(type,'hg')
     if size(terms,1)==1 && size(terms,2)==1                                  % Make terms request matrix

        n=ceil(sqrt(terms));

        tmat=ones(n,n);

-       if n^2>terms

-           tmat(end,end-ceil((n^2-terms)/2)+1:end)=0;

-           tmat(end-floor((n^2-terms)/2):end-1,end)=0;       

-       end

+       %if n^2>terms

+       %    tmat(end,end-ceil((n^2-terms)/2)+1:end)=0;

+       %    tmat(end-floor((n^2-terms)/2):end-1,end)=0;       

+       %end

    else

        tmat=terms;

    end

@@ -123,7 +134,7 @@ if strcmpi(type,'hg')
             m=t-1;

             if tmat(s,t)==1 

                 z2=HermiteGaussianE([l,m,q,lambda],domain(:,:,1),domain(:,:,2));

-                coeffs(s,t)=overlap(z1,conj(z2),domain,1,accuracy);

+                coeffs(s,t)=overlap(z1,conj(z2),domain,metric,accuracy);

             else

                 coeffs(s,t)=0;

             end

@@ -138,14 +149,14 @@ elseif strcmpi(type,'lg')
     if size(terms,1)==1 && size(terms,2)==1                                  % Make terms request matrix

        n=ceil( (1+sqrt(1+8*terms))/4 );

        tmat=ones(n,n,2); tmat(:,1,2)=0;

-       ndiff=2*n^2-n-terms;

-       if ndiff>0

-           rdiff=ceil(ndiff/2); ldiff=floor(ndiff/2);

-           tmat(end,end-ceil(rdiff/2)+1:end,1)=0;

-           tmat(end-floor(rdiff/2):end-1,end,1)=0;       

-           tmat(end,end-ceil(ldiff/2)+1:end,2)=0;

-           tmat(end-floor(ldiff/2):end-1,end,2)=0;  

-       end

+       %ndiff=2*n^2-n-terms;

+       %if ndiff>0

+       %    rdiff=ceil(ndiff/2); ldiff=floor(ndiff/2);

+       %    tmat(end,end-ceil(rdiff/2)+1:end,1)=0;

+       %    tmat(end-floor(rdiff/2):end-1,end,1)=0;       

+       %    tmat(end,end-ceil(ldiff/2)+1:end,2)=0;

+       %    tmat(end-floor(ldiff/2):end-1,end,2)=0;  

+       %end

        %disp(' '); dispmat(tmat(:,:,1)); disp(' '); dispmat(tmat(:,:,2)); disp(' '); disp(num2str(sum(sum(sum(tmat)))));

    else

        tmat=terms;

@@ -158,14 +169,14 @@ elseif strcmpi(type,'lg')
         for t=1:size(tmat,2)

             m=t-1;

             if tmat(s,t,1)==1                                               % coeff requested

-                z2=LaguerreGaussianE([p,m,q,lambda],domain(:,:,1),domain(:,:,2),'pol');

-                coeffs(s,t,1)=overlap(z1,conj(z2),domain,domain(:,:,1),accuracy);

+                z2=LaguerreGaussianE([p,m,q,lambda],domain(:,:,1),domain(:,:,2),coordtype);

+                coeffs(s,t,1)=overlap(z1,conj(z2),domain,metric,accuracy);

             else

                 coeffs(s,t,1)=0;

             end

             if tmat(s,t,2)==1                                               % coeff requested 

-                z2=LaguerreGaussianE([p,-m,q,lambda],domain(:,:,1),domain(:,:,2),'pol');

-                coeffs(s,t,2)=overlap(z1,conj(z2),domain,domain(:,:,1),accuracy);

+                z2=LaguerreGaussianE([p,-m,q,lambda],domain(:,:,1),domain(:,:,2),coordtype);

+                coeffs(s,t,2)=overlap(z1,conj(z2),domain,metric,accuracy);

             else

                 coeffs(s,t,2)=0;

             end