function plotBetter(A,j)

maxX = size(A,1);
maxY = size(A,2);
maxZ = size(A,3);

R = 3;   % increase resolution
L = 1*R;   % decay length
d = 2*L; % include so many points
s = 0.0;
for ix=-d:d
    for iy=-d:d
        for iz=-d:d
            d2 = ix*ix+iy*iy+iz*iz;
            s = s + exp(-d2/(L*L));
        end
    end
end
s;

maxXX = maxX*R;
maxYY = maxY*R;
maxZZ = maxZ*R;
B = zeros(maxXX, maxYY, maxZZ);
for x=1:maxX
%    fprintf('x= %d\n',x);
    fprintf('.');
    for y=1:maxY
        for z=1:maxZ
            if (A(x,y,z)==1) 
                for ix=-d:d
                    for iy=-d:d
                        for iz=-d:d
                            d2 = ix*ix+iy*iy+iz*iz;
                            f = exp(-d2/(L*L)) / s;
                            xx = (x-1)*R+ix+1;
                            yy = (y-1)*R+iy+1;
                            zz = (z-1)*R+iz+1;
                            xx = mod(xx-1,maxXX)+1;
                            yy = mod(yy-1,maxYY)+1;
%                            zz = mod(zz-1,maxZZ)+1;
%                            B(xx,yy,zz) = B(xx,yy,zz)+f;
                            if ((zz>=1) && (zz<=maxZZ)) B(xx,yy,zz) = B(xx,yy,zz)+f; end 
                        end
                    end
                end
            end
        end
    end
end
fprintf('\n');

%A
%B
%pause

clf
[XX,YY,ZZ] = meshgrid(0:maxXX-1,0:maxYY-1,0:maxZZ-1);
v = (max(max(max(B)))+min(min(min(B))))/2.0;
p = patch(isosurface(XX,YY,ZZ,B,v));
isonormals(XX,YY,ZZ,B,p);
set(p,'FaceColor','red','EdgeColor','none');
alpha(0.5);
daspect([1 1 1])
view(3); 
axis tight
axis([0,maxXX-1,0,maxYY-1,0,maxZZ-1]);
grid on
camlight 
lighting gouraud

if nargin > 1
    saveas(gcf, ['DLA3D_',sprintf('%03d',j),'.png'])
end

return

