% Shows Sinogram of harmonic functions
N = 256;
kval = 3*2*pi/N;
kang = 20;
cs = cos(kang); sn = sin(kang);
imi = sqrt(-1);
im = zeros(N);
for i = 1:N
    for j = 1:N
        im(i,j) = exp(imi*kval*(i*cs + j*sn));
    end
end
figure(1);clf;
subplot(2,2,1); imagesc(real(im));title('Real part of plane wave');
subplot(2,2,2); imagesc(imag(im));title('Imaginary part of plane wave');

%%

ang = [0:179];

g = radon(im,ang);

figure(1);
subplot(2,2,3); imagesc(real(g));title('Real part of sinogram');
subplot(2,2,4); imagesc(imag(g));title('Imaginary part of Sinogram')


%% now run through angles and frequencies

for kk = 1:10
    for kang = 1:20:180
        kval = kk*2*pi/N;
        cs = cos(kang); sn = sin(kang);
        im = zeros(N);
        for i = 1:N
            for j = 1:N
                im(i,j) = exp(imi*kval*(i*cs + j*sn));
            end
        end
        figure(1);clf;
        subplot(2,2,1); imagesc(real(im));title(['k = ',num2str(kk),' \theta=',num2str(kang),' (Real)']);
        subplot(2,2,2); imagesc(imag(im));title('(Imag)');
        
        g = radon(im,ang);
        
        subplot(2,2,3); imagesc(real(g));title('Real part of sinogram');
        subplot(2,2,4); imagesc(imag(g));title('Imaginary part of Sinogram');
        pause(1);
    end
end