function [S,Savg,x,y,xavg,yavg] =  blockAvg_better_1(Nx,Ny,Nxavg,Nyavg)
% For example, Nx = Ny = 1000, Nxavg = Nyavg = 25;
  
Lx = 10;
Ly = 10;


x = linspace(1,Lx,Nx);
y = linspace(1,Ly,Ny);

S = zeros(Nx,Ny);
for i = 1:length(x)
    for j = 1:length(y)
        yp = y(j);
        xp = x(i);
        S(i,j) = 5*cos((xp+yp)*2*pi) + 2*sin(2*xp*pi) + 2*cos(2*xp*pi);
    end
end

Savg = zeros(Nx,Ny);
yavg = zeros(floor(Ny/Nyavg),1 );
xavg = zeros(floor(Nx/Nxavg),1 );
for i = Nxavg:length(x)
    for j = Nyavg:length(y)
        
        if (mod(j,Nyavg)==0 && mod(i,Nxavg)==0)
            idown = i - Nxavg + 1;
            jdown = j - Nyavg + 1;
            iup = i;
            jup = j;
            
            iavg = floor(i/Nxavg);
            javg = floor(j/Nyavg);
            
            Savg(iavg,javg) = sum(sum(S(idown:iup,jdown:jup)))/...
                (Nxavg*Nyavg);
            xavg(iavg) = 0.5*(x(iup) + x(idown));
            yavg(javg) = 0.5*(y(jup) + y(jdown));
        end
        
    end
end
