high_res_cube/matlab-grid-generator/reduce_dx.m

function [dxg,dxc,dxf,dxv]=reduce_dx(dx,nratio);
% [dxg,dxc,dxf,dxv]=reduce_dx(dx,nratio);
%
% Sum dx on fine grid -> dxg,dxc,dxf,dxv
% nratio is ratio of grid sizes (e.g. nratio ~ dxg/dx)

if nratio==1
 dxg=dx;
 dxf=(dx(:,1:end-1)+dx(:,2:end))/2;
 dxv=(dx([end 1:end],:)+dx([1:end 1],:))/2;
 dxc=(dxf([end 1:end],:)+dxf([1:end 1],:))/2;
 return
end

[nxf]=size(dx,2);

if nratio/2 ~= floor(nratio/2)
 nratio
 error('nratio must be multiple of 2 to be able to reduce grid');
end

kg=1:nratio:nxf;
kc=1+nratio/2:nratio:nxf;
jg=1:nratio:nxf-1;
jc=[nxf-nratio/2 1+nratio/2:nratio:nxf];
dxg=dx(jg,kg);
dxf=dx(jg,kc);
dxc=dx(jc,kc);
dxv=dx(jc,kg);
for n=2:nratio;
 jg=mod(jg,nxf-1)+1;
 jc=mod(jc,nxf-1)+1;
 dxg=dxg+dx(jg,kg);
 dxf=dxf+dx(jg,kc);
 dxc=dxc+dx(jc,kc);
 dxv=dxv+dx(jc,kg);
end

% Force more symmetry
dxf=(dxf+dxf(end:-1:1,:))/2;
dxf=(dxf+dxf(:,end:-1:1))/2;
dxg=(dxg+dxg(end:-1:1,:))/2;
dxg=(dxg+dxg(:,end:-1:1))/2;
dxc=(dxc+dxc(end:-1:1,:))/2;
dxc=(dxc+dxc(:,end:-1:1))/2;
dxv=(dxv+dxv(end:-1:1,:))/2;
dxv=(dxv+dxv(:,end:-1:1))/2;
1	function [dxg,dxc,dxf,dxv]=reduce_dx(dx,nratio);
2	% [dxg,dxc,dxf,dxv]=reduce_dx(dx,nratio);
3	%
4	% Sum dx on fine grid -> dxg,dxc,dxf,dxv
5	% nratio is ratio of grid sizes (e.g. nratio ~ dxg/dx)
6
7	if nratio==1
8	dxg=dx;
9	dxf=(dx(:,1:end-1)+dx(:,2:end))/2;
10	dxv=(dx([end 1:end],:)+dx([1:end 1],:))/2;
11	dxc=(dxf([end 1:end],:)+dxf([1:end 1],:))/2;
12	return
13	end
14
15	[nxf]=size(dx,2);
16
17	if nratio/2 ~= floor(nratio/2)
18	nratio
19	error('nratio must be multiple of 2 to be able to reduce grid');
20	end
21
22	kg=1:nratio:nxf;
23	kc=1+nratio/2:nratio:nxf;
24	jg=1:nratio:nxf-1;
25	jc=[nxf-nratio/2 1+nratio/2:nratio:nxf];
26	dxg=dx(jg,kg);
27	dxf=dx(jg,kc);
28	dxc=dx(jc,kc);
29	dxv=dx(jc,kg);
30	for n=2:nratio;
31	jg=mod(jg,nxf-1)+1;
32	jc=mod(jc,nxf-1)+1;
33	dxg=dxg+dx(jg,kg);
34	dxf=dxf+dx(jg,kc);
35	dxc=dxc+dx(jc,kc);
36	dxv=dxv+dx(jc,kg);
37	end
38
39	% Force more symmetry
40	dxf=(dxf+dxf(end:-1:1,:))/2;
41	dxf=(dxf+dxf(:,end:-1:1))/2;
42	dxg=(dxg+dxg(end:-1:1,:))/2;
43	dxg=(dxg+dxg(:,end:-1:1))/2;
44	dxc=(dxc+dxc(end:-1:1,:))/2;
45	dxc=(dxc+dxc(:,end:-1:1))/2;
46	dxv=(dxv+dxv(end:-1:1,:))/2;
47	dxv=(dxv+dxv(:,end:-1:1))/2;