high_res_cube/matlab-grid-generator/reduce_E.m

function [Ec,Ez,Ev]=reduce_E(E,nratio);
% [Ec,Ez,Ev]=reduce_E(E,nratio);
%
% Sum E on fine grid -> Ec,Ez,Ev
% nratio is ratio of grid sizes (e.g. nratio ~ sqrt(Ec/E))

[nxf]=size(E,1)+1;

if nratio==1
 Ec=E;
 Ev=(E(:,[end 1:end])+E(:,[1:end 1]))/2;
 Ez=(Ev([end 1:end],:)+Ev([1:end 1],:))/2;
 return
end

if nratio/2 ~= floor(nratio/2)
 nratio
 error('nratio must be multiple of 2 to be able to reduce gid');
end
nx=(nxf-1)/nratio+1;
if nx ~= floor(nx)
 nx
 error('nxf/nratio must be an integer');
end

Ec=zeros(nx-1,nx-1);
Ev=zeros(nx-1,nx);
Ez=zeros(nx,nx);
kg=(1:nratio:nxf-1)-1;
kc=([nxf-nratio/2 1+nratio/2:nratio:nxf])-1;
for m=1:nratio;
 kg=mod(kg,nxf-1)+1;
 kc=mod(kc,nxf-1)+1;
 jg=(1:nratio:nxf-1)-1;
 jc=([nxf-nratio/2 1+nratio/2:nratio:nxf])-1;
 for n=1:nratio;
  jg=mod(jg,nxf-1)+1;
  jc=mod(jc,nxf-1)+1;
  Ec=Ec+E(jg,kg);
  Ev=Ev+E(jg,kc);
  Ez=Ez+E(jc,kc);
 end
end

% Force more symmetry
Ec=(Ec+Ec(end:-1:1,:))/2;
Ec=(Ec+Ec(:,end:-1:1))/2;
Ez=(Ez+Ez(end:-1:1,:))/2;
Ez=(Ez+Ez(:,end:-1:1))/2;
Ev=(Ev+Ev(end:-1:1,:))/2;
Ev=(Ev+Ev(:,end:-1:1))/2;

% Adjust corner Ez to a hexagon
Ez([1 end],[1 end])=Ez([1 end],[1 end])*(3/4);
1	function [Ec,Ez,Ev]=reduce_E(E,nratio);
2	% [Ec,Ez,Ev]=reduce_E(E,nratio);
3	%
4	% Sum E on fine grid -> Ec,Ez,Ev
5	% nratio is ratio of grid sizes (e.g. nratio ~ sqrt(Ec/E))
6
7	[nxf]=size(E,1)+1;
8
9	if nratio==1
10	Ec=E;
11	Ev=(E(:,[end 1:end])+E(:,[1:end 1]))/2;
12	Ez=(Ev([end 1:end],:)+Ev([1:end 1],:))/2;
13	return
14	end
15
16	if nratio/2 ~= floor(nratio/2)
17	nratio
18	error('nratio must be multiple of 2 to be able to reduce gid');
19	end
20	nx=(nxf-1)/nratio+1;
21	if nx ~= floor(nx)
22	nx
23	error('nxf/nratio must be an integer');
24	end
25
26	Ec=zeros(nx-1,nx-1);
27	Ev=zeros(nx-1,nx);
28	Ez=zeros(nx,nx);
29	kg=(1:nratio:nxf-1)-1;
30	kc=([nxf-nratio/2 1+nratio/2:nratio:nxf])-1;
31	for m=1:nratio;
32	kg=mod(kg,nxf-1)+1;
33	kc=mod(kc,nxf-1)+1;
34	jg=(1:nratio:nxf-1)-1;
35	jc=([nxf-nratio/2 1+nratio/2:nratio:nxf])-1;
36	for n=1:nratio;
37	jg=mod(jg,nxf-1)+1;
38	jc=mod(jc,nxf-1)+1;
39	Ec=Ec+E(jg,kg);
40	Ev=Ev+E(jg,kc);
41	Ez=Ez+E(jc,kc);
42	end
43	end
44
45	% Force more symmetry
46	Ec=(Ec+Ec(end:-1:1,:))/2;
47	Ec=(Ec+Ec(:,end:-1:1))/2;
48	Ez=(Ez+Ez(end:-1:1,:))/2;
49	Ez=(Ez+Ez(:,end:-1:1))/2;
50	Ev=(Ev+Ev(end:-1:1,:))/2;
51	Ev=(Ev+Ev(:,end:-1:1))/2;
52
53	% Adjust corner Ez to a hexagon
54	Ez([1 end],[1 end])=Ez([1 end],[1 end])*(3/4);