anim=0;
seed=1;  %change this to get different random realization
ns=1;    %spectral index with ns=1 scale invariant
sigma=0; %smoothing scale

sz=512;
mk=zeros(sz,sz);


if (anim==1)
 loops=40;
 aviobj = avifile('z:\power_ns.avi','colormap',jet,'fps',2,'compression','Cinepak');
 winsize = get(gcf,'Position');
 winsize(1:2) = [0 0];
else
    loops=1;
end

for loop=1:loops

    randn('state',seed);
    
    if (anim==1)
    ns=-2+ 6*(loop-1)/loops;
    end

for kx=1:sz/2
 for ky=1:sz
      ky2=mod(sz-ky+1,sz) + 1;
      kky=min(ky,ky2);
      k = sqrt((kx-1)^2+(kky-1)^2);
      if (k>0)
           if (sigma==0)
               amp=k^((ns-1)/2)/k;
           else
               amp=k^((ns-1)/2)/k *exp(-((k*sigma))^2/2);
           end
      else
          amp=0;
      end;
    if (kx==1 && ky2==ky)
      mk(kx,ky) = amp*randn*sqrt(2); 
     else
      mk(kx,ky) = amp*complex(randn,randn);    
    end  
 end
end


for kx=sz/2+1:sz
 for ky=1:sz
      kx2=mod(sz-kx+1,sz) + 1;
      ky2=mod(sz-ky+1,sz) + 1;
      mk(kx,ky) =conj( mk(kx2,ky2));
 end
end

m=ifft2(mk,'symmetric');
m=128+ m*256/max(max(abs(m)));

image(m);
axis off          % Remove axis ticks and numbers
axis image 
colormap(jet);

if (anim==1)
aviobj = addframe(aviobj,m);
end

end %loops

if (anim==1)
 aviobj = close(aviobj);
else
    cols=colormap;
%imwrite(m,cols,'z:\scale-inv.bmp','bmp');
end