mat = 64;
m2 = mat/2;
[x y] = meshgrid(-m2:m2-1);
r = sqrt(x.^2+y.^2);
phantom = r<0.5*m2;
data = ft(phantom);

rpow = r.^(0.8); % assume time is proportional to sqrt(r)
rpowmax = rpow(1,m2);
rpow = (rpow<rpowmax) .* (rpow / rpowmax);

phmax = 2*pi; % max phase rotation across data acquisition
    
figure(777);
subplot(2,3,4);
imagesc(phantom); axis off; axis image;
subplot(2,3,1);
imagesc(real(data)); axis off; axis image;
kmod = exp(j*x*phmax/mat);
subplot(2,3,2);
imagesc(real(kmod)); axis off; axis image;
subplot(2,3,5);
imagesc(abs(ift(kmod))); axis off; axis image;
subplot(2,3,3);
ph = ift((kmod.*data));
plot(phantom(m2,:),'b');
hold on;
plot(real(ph(m2,:)),'r'); axis([0 mat+1 -0.1 1.1]);
%imagesc(real(kmod.*data)); axis off;
subplot(2,3,6);
imagesc(abs(ph)); axis off; axis image;

figure(888);
subplot(2,3,4);
imagesc(phantom); axis off; axis image;
subplot(2,3,1);
imagesc(real(data)); axis off; axis image;
kmod = exp(j*rpow*phmax);
subplot(2,3,2);
imagesc(real(kmod)); axis off; axis image;
subplot(2,3,5);
imagesc(abs(ift(kmod))); axis off; axis image;
subplot(2,3,3);
ph = ift((kmod.*data));
plot(phantom(m2,:),'b');
hold on;
plot(real(ph(m2,:)),'r'); axis([0 mat+1 -0.1 1.1]);
%imagesc(real(kmod.*data)); axis off;
subplot(2,3,6);
imagesc(abs(ph)); axis off; axis image;
