%Generates results and plots for Chapter 4 Tomography example 4.1
clear

t=sqrt(2);

%load the system matrix for these ray paths

G = [1,0,0,1,0,0,1,0,0;
     0,1,0,0,1,0,0,1,0;
     0,0,1,0,0,1,0,0,1;
     1,1,1,0,0,0,0,0,0;
     0,0,0,1,1,1,0,0,0;
     0,0,0,0,0,0,1,1,1;
     t,0,0,0,t,0,0,0,t;
     0,0,0,0,0,0,0,0,t];

%get the singular values
[U,s,V] = svd(G);

disp('Singular values:')
diag(s)

disp('System rank:')

[m,n]=size(G);

r=rank(G);

disp('System null space vectors')
V(:,r+1:n)


disp('System null space vectors reshaped into matrices')
m01=reshape(V(:,r+1),3,3)
m02=reshape(V(:,r+2),3,3)

figure(1)
bookfonts;
colormap('gray')
imagesc(m01)
H=colorbar;
set(H,'Fontsize',18);
set(gca,'xtick',[1,2,3]);
set(gca,'ytick',[1,2,3]);
xlabel('j')
ylabel('i')
%print -deps2 fig4_1_1.eps

figure(2)
bookfonts;
colormap('gray')
imagesc(m02)
H=colorbar;
set(H,'Fontsize',18);
set(gca,'xtick',[1,2,3]);
set(gca,'ytick',[1,2,3]);
xlabel('j')
ylabel('i')
%print -deps2 fig4_1_2.eps

disp('System data space null vector')
U(:,m)

%model resolution matrix

Vp=V(:,1:r);
Rm=Vp*Vp'

CLIM=[min(min(Rm)),max(max(Rm))];

figure(3)
bookfonts;
colormap('gray')
imagesc(Rm,CLIM)
H=colorbar;
set(H,'Fontsize',18);
set(gca,'xtick',[1,2,3,4,5,6,7,8,9]);
set(gca,'ytick',[1,2,3,4,5,6,7,8,9]);
xlabel('j')
ylabel('i')
%print -deps2 fig4_1_3.eps

Rmdiag=reshape(diag(Rm),3,3);

figure(4)
colormap('gray')
bookfonts;
imagesc(Rmdiag,CLIM)
H=colorbar;
set(H,'Fontsize',18);
set(gca,'xtick',[1,2,3]);
set(gca,'ytick',[1,2,3]);
xlabel('j')
ylabel('i')
%print -deps2 fig4_1_4.eps

%spike resolution test

mtest=zeros(n,1);
mtest(5)=1;

dtest=G*mtest;

disp('model spike and predicted data')

mtest

dtest

disp('recovered model from the spike test')

mdagger=pinv(G)*dtest;

CLIM=[min(min([mtest;mdagger])),max(max([mtest;mdagger]))];

figure(5)
colormap('gray')
bookfonts
imagesc(reshape(mtest,3,3),CLIM);
H=colorbar;
set(H,'Fontsize',18);
set(gca,'xtick',[1,2,3]);
set(gca,'ytick',[1,2,3]);
xlabel('j')
ylabel('i')
%print -deps2 fig4_1_5.eps

figure(6)
colormap('gray')
bookfonts
imagesc(reshape(mdagger,3,3),CLIM);
H=colorbar;
set(H,'Fontsize',18);
set(gca,'xtick',[1,2,3]);
set(gca,'ytick',[1,2,3]);
xlabel('j')
ylabel('i')
%print -deps2 fig4_1_6.eps





