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);