osse/utils/cyl2cartuv.m

function [u,v] = cyl2cartuv(thetav,rhov,xi,yi,varargin)
% [u,v]=cyl2cartuv(thetav,rhov,xi,yi);
%
% Re-grids model output in cylindrical coords to cartesian.
%  c     is a 2-D or 3-D scalar or z-vector field
%  xi,yi are vectors of the new regular lat-lon grid to interpolate to.
%  z     is the interpolated data with dimensions of size(xi) by size(yi).
%  theta=0 is at 12 o'clock.
%
% e.g.
% >> t=rdmds('Ttave.0000513360');
% >> xi=-179:2:180;yi=-89:2:90;
% >> ti=cyl2cart(x,y,t,xi,yi);
%
% $Header: /u/gcmpack/MITgcm_contrib/osse/utils/cyl2cartuv.m,v 1.1 2004/06/08 14:14:44 afe Exp $

if ~isequal(size(thetav),size(rhov))
  error('Theta and rho vector arrays must be same size');
end

%work out mappings of polar to cartesian
NN=size(thetav);
[ntheta nrho nz]=size(thetav);
[RHO,THETA,NZ] = meshgrid(1:nrho,-pi+2*pi/ntheta:2*pi/ntheta:pi,1:nz);
[x,y] = pol2cart(THETA(:,:,1),RHO(:,:,1));
%[nx ny nz]=size(c);
nx=ntheta;ny=nrho;

% break out components
%vv=-(thetav.*-cos(THETA)+rhov.*sin(THETA));
%uv=-(thetav.*-sin(THETA)+rhov.*cos(THETA));
uv=thetav.*cos(THETA)+rhov.*sin(THETA);
vv=thetav.*-sin(THETA)+rhov.*cos(THETA);
%uv=thetav.*cos(THETA); %+rhov.*sin(THETA);
%vv=thetav.*-sin(THETA); %+rhov.*cos(THETA);
%uv=rhov.*sin(THETA);
%vv=rhov.*cos(THETA);


X=reshape(x,[1 nx*ny]);
Y=reshape(y,[1 nx*ny]);
del=griddata_preprocess(Y,X,yi,xi',varargin{:});

for k=1:nz;
 UV=reshape(uv(:,:,k),[1 nx*ny]);
 VV=reshape(vv(:,:,k),[1 nx*ny]);
u(:,:,k)=griddata(Y,X,UV,yi,xi',varargin{:});
v(:,:,k)=griddata(Y,X,VV,yi,xi',varargin{:});
% z(:,:,k)=griddata_fast(del,[C C(il) C(ig)],varargin{:});
end % k

% Split vertical and time dimensions
if size(NN,2)>2
u=reshape(u,[size(u,1) size(u,2) NN(3:end)]);
v=reshape(v,[size(v,1) size(v,2) NN(3:end)]);
end
1	afe	1.1	function [u,v] = cyl2cartuv(thetav,rhov,xi,yi,varargin)
2			% [u,v]=cyl2cartuv(thetav,rhov,xi,yi);
3			%
4			% Re-grids model output in cylindrical coords to cartesian.
5			% c is a 2-D or 3-D scalar or z-vector field
6			% xi,yi are vectors of the new regular lat-lon grid to interpolate to.
7			% z is the interpolated data with dimensions of size(xi) by size(yi).
8	afe	1.2	% theta=0 is at 12 o'clock.
9	afe	1.1	%
10			% e.g.
11			% >> t=rdmds('Ttave.0000513360');
12			% >> xi=-179:2:180;yi=-89:2:90;
13			% >> ti=cyl2cart(x,y,t,xi,yi);
14			%
15	afe	1.2	% $Header: /u/gcmpack/MITgcm_contrib/osse/utils/cyl2cartuv.m,v 1.1 2004/06/08 14:14:44 afe Exp $
16	afe	1.1
17			if ~isequal(size(thetav),size(rhov))
18			error('Theta and rho vector arrays must be same size');
19			end
20
21			%work out mappings of polar to cartesian
22			NN=size(thetav);
23			[ntheta nrho nz]=size(thetav);
24			[RHO,THETA,NZ] = meshgrid(1:nrho,-pi+2pi/ntheta:2pi/ntheta:pi,1:nz);
25			[x,y] = pol2cart(THETA(:,:,1),RHO(:,:,1));
26			%[nx ny nz]=size(c);
27			nx=ntheta;ny=nrho;
28
29			% break out components
30			%vv=-(thetav.-cos(THETA)+rhov.sin(THETA));
31			%uv=-(thetav.-sin(THETA)+rhov.cos(THETA));
32			uv=thetav.cos(THETA)+rhov.sin(THETA);
33			vv=thetav.-sin(THETA)+rhov.cos(THETA);
34			%uv=thetav.cos(THETA); %+rhov.sin(THETA);
35			%vv=thetav.-sin(THETA); %+rhov.cos(THETA);
36			%uv=rhov.*sin(THETA);
37			%vv=rhov.*cos(THETA);
38
39
40			X=reshape(x,[1 nx*ny]);
41			Y=reshape(y,[1 nx*ny]);
42			del=griddata_preprocess(Y,X,yi,xi',varargin{:});
43
44			for k=1:nz;
45			UV=reshape(uv(:,:,k),[1 nx*ny]);
46			VV=reshape(vv(:,:,k),[1 nx*ny]);
47			u(:,:,k)=griddata(Y,X,UV,yi,xi',varargin{:});
48			v(:,:,k)=griddata(Y,X,VV,yi,xi',varargin{:});
49			% z(:,:,k)=griddata_fast(del,[C C(il) C(ig)],varargin{:});
50			end % k
51
52			% Split vertical and time dimensions
53			if size(NN,2)>2
54			u=reshape(u,[size(u,1) size(u,2) NN(3:end)]);
55			v=reshape(v,[size(v,1) size(v,2) NN(3:end)]);
56			end