high_res_cube/matlab-grid-generator/calc_fvgrid.m

function [dx,dy,E] = calc_fvgrid(lx,ly)
% Calculates finite volume grid info (dx,dy,A) for conformal cubic grid
% with curvilinear coordinates (lx,ly)
%
% Meant to be used for single quadrant of tile but does work for full tile

nxf=size(lx,1);

%j=1:nx-1;jm1=j-1;jm1(1)=nx-1;
%J=1:nx;Jm1=J-1;Jm1(1)=nx-1;
%J=j;

% 3D coordinates on unit sphere
[lX,lY,lZ]=map_xy2xyz(lx,ly);

% Symmetrize
%lX=(lX-lX(end:-1:1,:))/2; lX=(lX+lX(:,end:-1:1))/2;
%lY=(lY+lY(end:-1:1,:))/2; lY=(lY-lY(:,end:-1:1))/2;
%lZ=(lZ+lZ(end:-1:1,:))/2; lZ=(lZ+lZ(:,end:-1:1))/2;

% Symmetry?
%disp('calc_grid: symmetry tests of 3D coordinates');
%stats( lX+lX(end:-1:1,:) )
%stats( lX-lX(:,end:-1:1) )
%stats( lY-lY(end:-1:1,:) )
%stats( lY+lY(:,end:-1:1) )
%stats( lZ-lZ(end:-1:1,:) )
%stats( lZ-lZ(:,end:-1:1) )
%stats( lZ-lZ' )

% Use "vector" data type
Vertices=coord2vector(lX,lY,lZ);
dx=angle_between_vectors(Vertices(1:end-1,:,:),Vertices(2:end,:,:));
dy=angle_between_vectors(Vertices(:,1:end-1,:),Vertices(:,2:end,:));
E=excess_of_quad(Vertices(1:end-1,1:end-1,:),Vertices(2:end,1:end-1,:), ...
                 Vertices(2:end,2:end,:),Vertices(1:end-1,2:end,:));

% Force some symmetry (probably does nothing)
dx=(dx+dy')/2;
dy=dx';
E=(E+E')/2;

% Use symmetry to fill second octant
for j=2:nxf;
 dx(1:j-1,j)=dy(j,1:j-1)';
end
for j=1:nxf-1;
 dy(1:j,j)=dx(j,1:j)';
end
for j=2:nxf-1;
 E(1:j-1,j)=E(j,1:j-1)';
end
1	function [dx,dy,E] = calc_fvgrid(lx,ly)
2	% Calculates finite volume grid info (dx,dy,A) for conformal cubic grid
3	% with curvilinear coordinates (lx,ly)
4	%
5	% Meant to be used for single quadrant of tile but does work for full tile
6
7	nxf=size(lx,1);
8
9	%j=1:nx-1;jm1=j-1;jm1(1)=nx-1;
10	%J=1:nx;Jm1=J-1;Jm1(1)=nx-1;
11	%J=j;
12
13	% 3D coordinates on unit sphere
14	[lX,lY,lZ]=map_xy2xyz(lx,ly);
15
16	% Symmetrize
17	%lX=(lX-lX(end:-1:1,:))/2; lX=(lX+lX(:,end:-1:1))/2;
18	%lY=(lY+lY(end:-1:1,:))/2; lY=(lY-lY(:,end:-1:1))/2;
19	%lZ=(lZ+lZ(end:-1:1,:))/2; lZ=(lZ+lZ(:,end:-1:1))/2;
20
21	% Symmetry?
22	%disp('calc_grid: symmetry tests of 3D coordinates');
23	%stats( lX+lX(end:-1:1,:) )
24	%stats( lX-lX(:,end:-1:1) )
25	%stats( lY-lY(end:-1:1,:) )
26	%stats( lY+lY(:,end:-1:1) )
27	%stats( lZ-lZ(end:-1:1,:) )
28	%stats( lZ-lZ(:,end:-1:1) )
29	%stats( lZ-lZ' )
30
31	% Use "vector" data type
32	Vertices=coord2vector(lX,lY,lZ);
33	dx=angle_between_vectors(Vertices(1:end-1,:,:),Vertices(2:end,:,:));
34	dy=angle_between_vectors(Vertices(:,1:end-1,:),Vertices(:,2:end,:));
35	E=excess_of_quad(Vertices(1:end-1,1:end-1,:),Vertices(2:end,1:end-1,:), ...
36	Vertices(2:end,2:end,:),Vertices(1:end-1,2:end,:));
37
38	% Force some symmetry (probably does nothing)
39	dx=(dx+dy')/2;
40	dy=dx';
41	E=(E+E')/2;
42
43	% Use symmetry to fill second octant
44	for j=2:nxf;
45	dx(1:j-1,j)=dy(j,1:j-1)';
46	end
47	for j=1:nxf-1;
48	dy(1:j,j)=dx(j,1:j)';
49	end
50	for j=2:nxf-1;
51	E(1:j-1,j)=E(j,1:j-1)';
52	end