function varargout=localization2D(XY,N,J,method,Nqx,Nqy,dxdy,XYP) % [G,H,V,K,XYP,XY,A,c11,cmn]=LOCALIZATION2D(XY,N,J,method,Nqx,Nqy,dxdy,XYP) % % Calculates spatiospectrally concentrated Slepian functions in the % two-dimensional Cartesian plane at the desired spatial points, for an % abitrarily shaped region. These are windows to be applied to certain % data of which you know the specific coordinates or the x and y % spacings. If you don't, put in a reasonable sampling interval in km. % % INPUT: % % XY [X(:) Y(:)] coordinates of a closed curve [default: circle] % N Shannon number [default: 20] % J Number of eigentapers [default: 2*N] % method 'GL' by Gauss-Legendre integration [default] % 'RS' by equal-area Riemann summation % Nqx Number of x quadrature points [dx-based for 'RS' / 32 for 'GL'] % Nqy Number of y quadrature points [dy-based for 'RS' / 32 for 'GL'] % dxdy [dx dy] the x and y spacing for the return grid OR % XYP [X(:) Y(:)] evaluation points [default: rectangle inscribing region] % c11,cmn Coordinates of the grid corner points for plotting % % OUTPUT: % % G The first J space-concentrated eigentapers at XYP [size(XYP) J] % H The first J spacelimited eigentapers at XYP [size(XYP) J] % V The eigenvalues of the tapers [vector] % K The circular bandlimit wavenumber for this problem % XYP Return coordinates for H and G [unwrapped two-column matrix] % XY The spatial curve - again % A The approximate area of the domain enclosed by the points % % EXAMPLES: % % localization2D('demo1'); % For the circle % localization2D('demo2'); % For a random blob % localization2D('demo3'); % For a profile through the circle % localization2D('demo4'); % For a profile through the blob % localization2D('demo1',[],[],'RS'); % localization2D('demo1',[],[],'GL'); % localization2D('demo2',[],[],'RS'); % localization2D('demo2',[],[],'GL'); % localization2D('demo4',[],[],'RS'); % localization2D('demo4',[],[],'GL'); % localization2D('demo5') % Compare eigenvalues for circle % localization2D('demo6') % Compare eigenvalues for blob % localization2D('demo7') % For a square % % SEE ALSO: % % SWDISK, KERNELC2D, LOCALIZATION, KERNELC, DEVILLIERS % % Last modified by dongwang-at-princeton.edu, 05/14/2010 % Last modified by fjsimons-at-alum.mit.edu, 06/17/2010 % Supply defaults - make sure to hit the diameter exactly defval('circn',41) % Note that the discretization of the boundary curve makes for the % integration accuracy of the circle in the Gauss-Legendre case, and not % so much the number of nodes defval('XY',[cos(linspace(0,2*pi,circn)) ; sin(linspace(0,2*pi,circn))]') if ~isstr(XY) defval('N',20) defval('method','GL'); % Make sure the XY of the curve has two columns if ~isempty(find(size(XY,2)==2)) if size(XY,1)length(QinR) error('Increase quadrature resolution if you want this many eigenvalues') end % Report on the area using these approximations A=sum(w(:)); disp(sprintf('Approximate area is %8.3f',A)) % The wavenumber of the circular bandlimitation K=sqrt(4*pi*N/A); disp(sprintf('Starting kernel calculation of size %i x %i',... length(QinR),length(QinR))) % Now calculate the integral kernel at the same unwrapped set of % quadrature points in both of its slots D=kernelc2D([QX(QinR) QY(QinR)],[QX(QinR) QY(QinR)],K); % Square root of weight matrix - does precomputing it actually save time? sqrtW=sqrt(diag(w(:))); % Calculate the eigenvectors, [i.e. G and H on the q-grid] and eigenvalues % Here I'd substitute EIGS in case you don't want the whole set t0=clock; if J