% simple tomography
clear all;

% want to know m on Nx by Ny grid
Nx = 24;
xmin = -24.0;
xmax =  24.0;
dx=(xmax-xmin)/(Nx-1);
x = xmin+dx*[0:Nx-1]';
Dx = xmax-xmin;

Ny = 24;
ymin = -24.0;
ymax =  24.0;
dy=(ymax-ymin)/(Ny-1);
y = ymin+dy*[0:Ny-1]';
Dy = ymax-ymin;

% build true model
S = ones(Nx,Ny);
S( round(Nx/4) : round(3*Nx/4), round(Ny/4) : round(3*Ny/4) )=2.0;
S( round(Nx/8) : round(Nx/4) , round(Ny/8) : round(Ny/4) )=3.0;
S( round(3*Nx/8) : round(5*Nx/8), round(3*Ny/8) : round(5*Ny/8) )=4.0;
S( round(1*Nx/8) : round(2*Nx/8), round(5*Ny/8) : round(7*Ny/8) )=1.5;
S( round(5*Nx/8) : round(7*Nx/8), round(5*Ny/8) : round(7*Ny/8) )=5.0;
S( round(7*Nx/16) : round(9*Nx/16), round(7*Ny/16) : round(9*Ny/16) )=1.0;

% rearrage model into vector
% order or model parameters in vector: index = (ix-1)*Ny + iy
mtrue = zeros(Nx*Ny,1);
rowindex = zeros(Nx*Ny,1);
colindex = zeros(Nx*Ny,1);
for ix = [1:Nx]
for iy = [1:Ny]
    q = (ix-1)*Ny + iy;
    mtrue(q)=S(ix,iy);
    rowindex(q)=ix;
    colindex(q)=iy;
end
end
    
% plot true model
figure(1);
clf;
set(gca, 'YDir', 'reverse' );
imagesc( [xmin, xmax], [ymin, ymax], S);
hold on;
colorbar;
xlabel('y');
ylabel('x');
title( 'true model' );

% sources and receivers
Nside = 9;
Ns = 4*Nside;
xs = zeros(Ns, 2);
% left side
xs(1:Nside,1)=xmin+dx/10;
xs(1:Nside,2)=ymin+Dy*[1:Nside]'/(Nside+1);
% right side
xs(Nside+1:2*Nside,1)=xmax-dx/10;
xs(Nside+1:2*Nside,2)=ymin+Dy*[1:Nside]'/(Nside+1);
% bottom side
xs(2*Nside+1:3*Nside,1)=xmin+Dx*[1:Nside]'/(Nside+1);
xs(2*Nside+1:3*Nside,2)=ymin+dy/10;
% bottom side
xs(3*Nside+1:4*Nside,1)=xmin+Dx*[1:Nside]'/(Nside+1);
xs(3*Nside+1:4*Nside,2)=ymax-dy/10;

% plot receivers
plot( xs(:,2), xs(:,1), 'w+' );

% source receiver pairs
N = Ns*(Ns-1)/2;
ray = zeros(N,4);
k=1;
for i = [1:Ns]
    for j = [i+1:Ns]
        ray(k,1) = xs(i,1);
        ray(k,2) = xs(i,2);
        ray(k,3) = xs(j,1);
        ray(k,4) = xs(j,2);
        k=k+1;
    end
end

% plot true model and rays
figure(2);
clf;
set(gca, 'YDir', 'reverse' );
imagesc( [xmin, xmax], [ymin, ymax], S);
hold on;
colorbar;
xlabel('y');
ylabel('x');
title( 'true model with rays' );
plot( xs(:,2), xs(:,1), 'w+' );
for k = [1:N]
    plot( [ray(k,2) ray(k,4)], [ray(k,1), ray(k,3)], 'w' );
    end
    
% build data kernel
% order or model parameters in vector: index = (ix-1)*Ny + iy
M = Nx*Ny;
G=zeros(N,M);
Nr = 2000;
for k = [1:N]
    x1 = ray(k,1);
    y1 = ray(k,2);
    x2 = ray(k,3);
    y2 = ray(k,4);
    r = sqrt( (x2-x1)^2 + (y2-y1)^2 );
    dr = r/Nr;
% I use a sloppy way of computing the length of the ray
% in each pixel.  I subdivide the ray into Nr pieces, and
% assign each piece to exactly one pixel, the one its
% closest to
    for i = [1:Nr]
        x = x1 + (x2-x1)*i/Nr;
        y = y1 + (y2-y1)*i/Nr;
        ix = 1+floor( (Nx-1)*(x-xmin)/Dx );
        iy = 1+floor( (Ny-1)*(y-ymin)/Dy );
        q = (ix-1)*Ny + iy;
        G(k,q) = G(k,q) + dr;
    end
end

% compute true travel times and plot them
d = G*mtrue;

% plot true travel times on a 2D diagram
% arranged by the rays's closest approach to the origin
% and its orientation

% want to plot traveltimes on NR by Ntheta grid
NR = 100;
Rmin = 0.0;
Rmax = sqrt((Dx/2)^2+(Dy/2)^2);
dR = (Rmax-Rmin)/(NR-1);

Ntheta = 100;
thetamin = -180.0;
thetamax = 180.0;
dtheta = (thetamax-thetamin)/(Ntheta-1);

% build traveltime grid for plotting purposes
TT = zeros(NR,Ntheta);
s = zeros(N,1);
R = zeros(N,1);
theta = zeros(N,1);
for k = [1:N]
    x1 = ray(k,1);
    y1 = ray(k,2);
    x2 = ray(k,3);
    y2 = ray(k,4);
    % compute closest approach to origin and angle
    xd = (x2-x1);
    yd = (y2-y1);
    s(k) = -(x1*xd+y1*yd)/(xd^2+yd^2);
    xs = x1+xd*s(k);
    ys = y1+yd*s(k);
    R(k) = sqrt( xs^2 + ys^2 );
    theta(k) = (180/pi)*atan2(ys,xs);
    iR = 1+round((R(k)-Rmin)/dR);
    if( iR<1 )
        iR=1;
    elseif ( iR>NR )
        ir=NR;
    end
    itheta = 1+round((theta(k)-thetamin)/dtheta);
    if( itheta<1 )
        itheta=1;
    elseif ( itheta>Ntheta )
        itheta=Ntheta;
    end
    TT(iR,itheta)=d(k);
    % out in data for reverse orientation of ray, for symmetery
    theta2 = (180/pi)*atan2(-ys,-xs);
    itheta = 1+round((theta2-thetamin)/dtheta);
    if( itheta<1 )
        itheta=1;
    elseif ( itheta>Ntheta )
        itheta=Ntheta;
    end
    TT(iR,itheta)=d(k);

end

figure(3);
clf;
set(gca, 'YDir', 'reverse' );
imagesc( [thetamin, thetamax], [Rmin, Rmax], TT);
hold on;
colorbar;
xlabel('theta');
ylabel('R');
title( 'traveltimes' );

% inversion of data by damped least squares
mest=zeros(M,1);
epsilon = 1.0e-4;
GTG = G'*G;
for i = [1:M]
    GTG(i,i) = GTG(i,i)+epsilon;
end
GTd = (G'*d);
mest = GTG\GTd;

% rearrange m back into image
Sest = zeros(Nx,Ny);
for k = [1:M]
    ix=rowindex(k);
    iy=colindex(k);
    Sest(ix,iy)=mest(k);
end

% plot estimated image
% plot true model
figure(4);
clf;
set(gca, 'YDir', 'reverse' );
imagesc( [xmin, xmax], [ymin, ymax], Sest);
hold on;
colorbar;
xlabel('y');
ylabel('x');
title( 'estimated model' );