%1) create a mesh
rad = 25;   % mesh radius [mm]
nsect = 6;  % number of sectors
nring = 32; % number of rings
nbnd = 2;   % number of boundary rings

[vtx,idx,eltp] = mkcircle(rad,nsect,nring,nbnd);
% create the mesh geometry
mesh = toastMesh (vtx,idx,eltp);
% create the mesh object

[vtx,idx,eltp] = mesh.Data;
nnode = mesh.NodeCount;
nel = mesh.ElementCount;

%visualize the mesh
figure(1)
mesh.Display

%2) optical coeffs
%a)homogeneous distribution
mua_bkg = 0.01;
mus_bkg = 1.0;
ref_bkg = 1.4;
nnd = mesh.NodeCount;
mua = ones(nnd,1) * mua_bkg;
mus = ones(nnd,1) * mus_bkg;
ref = ones(nnd,1) * ref_bkg;

%b)source-detector locations
nq = 16;
for i=1:nq
   phi_q = 2*pi*(i-1)/nq;
   Q(i,:) = rad * [cos(phi_q) sin(phi_q)];
   phi_m = 2*pi*(i-0.5)/nq;
   M(i,:) = rad * [cos(phi_m) sin(phi_m)];
end
mesh.SetQM(Q,M); %attach the source-detector locations to the mesh

%visualization
hold on
plot(Q(:,1),Q(:,2),'ro','MarkerFaceColor','r');
plot(M(:,1),M(:,2),'bx','MarkerFaceColor','b');

%c)source and boundary projection vectors
qvec = mesh.Qvec('Neumann','Gaussian',2); %incoming flux
mvec = mesh.Mvec('Gaussian',2); %boundary profiles

%d)run the forward solver
K = dotSysmat(mesh,mua,mus,ref,0);
Phi = Kqvec;
Y = mvec.' * Phi;

#For larger problems, we have to switch to an iterative scheme such as bicgstab or gmres.

figure(2)
imagesc(log(Y));
xlabel('source index q');
ylabel('detector index m');
axis equal tight;
colorbar

figure(3)
hold on
angle = 360/32:360/16:360;
for i=1:size(Y,2)
   ywrap = [Y(i:end,i); Y(1:i-1,i)];
   plot(angle,log(ywrap),'o-');
end
axis([0 360 -13 -2]);
xlabel('angular source-detector separation');
ylabel('log intensity');