|
%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');
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-10-20 11:18
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社