Plot refletance

[Giovanna Tramontin]

Hello professor Fang, 

I am running a Monte Carlo to obtain the reflectance in a mesh slab geometry using 

 cfg.issaveref = 1;

cfg.isreflect = 1;

 cfg.outputtype = 'flux';

 [output] = mmclab(cfg);

I get as an output "output.dref", which is the reflectance data of each face of the surface of my mesh. And I retrieve the faces in the surface of the mesh using: 

 faces = faceneighbors(cfg.elem, 'rowmajor');

Also, I take the cw data: 

data = output.dref;
cwdata = squeeze(sum(data, ndims(data)));  

But I do not know how to plot this. 

 "cwdata" has the dimensions 4380x1, "faces" has the dimensions 4380x3 and is in terms of the indices of the nodes that connects each face. 

I would like to see the reflectance for the surface where z =10 but I do not know how to plot this and I did not find any similar example. Can you help me? Following below this email is a short version of the code. 

Best wishes,

Giovanna

%

clear path;
clear cfg;

%% Geometry:

% Slab properties
s0 = [0 0 0]; % Initial point of the slab
sf = [50 50 10]; % End point of the slab
p0 = [1 1 1];   % A point inside the slab
node_size = 2.0;
nd =  2.0;   %number of divisions

[node, face, elem] = meshabox(s0, sf, node_size, nd);
[cfg.node, cfg.elem] = surf2mesh(node, face, s0,sf, 1, nd, p0, [], 0);

cfg.elemprop=cfg.elem(:,5);
cfg.elem=cfg.elem(:,1:4);

g1 = 0;
n1 = 1.37;
mua1 = 0.003;  %mm^-1
mus1 = 1;    %mm^-1 - put here the value of mus'

cfg.prop = [0, 0, 1, 1; mua1, mus1, g1, n1];

cfg.unitinmm = 1;
cfg.basisorder = 1;
cfg.isreflect = 1; %Whether to do reflection/transmission; 0 no, 1 yes.
cfg.nphoton = 1e8;
cfg.gpuid  = 1;

cfg.tsart = 0;                
cfg.tend = 10e-9;            
cfg.tstep = 5e-10; %use 1024 bins, increase bins = lose resolution            

cfg.debuglevel = 'TP';

 cfg.issaveref = 1;

%% Source

%Source position (xyz)
cfg.srcpos = [0 0 10];  % one corner of the illumination area
%Source direction
cfg.srcdir = [0 0 -1];
pat.patternsize = 50.0;
pat.ky=0;                       % wave number in the y-dir
pat.yphase=0;  

%% Create light patterns and run simulation

for kx = 3 %[0: 1: 24]
    for xphase = [0, 2*pi/3, 4*pi/3]

        kx = kx;
        xphase = xphase;

        % Create light patterns
       
        method='elem';
        cfg.srctype='fourier';      % define an SFDI source
cfg.srcparam1=[patternsize 0 0 kx+xphase/(2*pi)];   % kx is k-number in x direction
cfg.srcparam2=[0 patternsize 0 ky+yphase/(2*pi)];
       
        cfg.outputtype = 'flux';
        [output] = mmclab(cfg);

        %% Plot data

        data = output.dref;
        cwdata = squeeze(sum(data, ndims(data)));  

        faces = faceneighbors(cfg.elem, 'rowmajor');
       
        surfaceelem =edgeneighbors(faces); %elements around this faces
       
    end
end

[Qianqian Fang]

hi Giovanna,

please use the below sample code provided in the basic mmclab example

https://github.com/fangq/mmc/blob/v2024.2/mmclab/example/demo_mmclab_basic.m#L37-L43

the tick is to set 'cdata' of the triangles in a patch object.

Qianqian

--
You received this message because you are subscribed to the Google Groups "mmc-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mmc-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mmc-users/27d3d4ce-e53b-4f8a-a365-5a8de46dc492n%40googlegroups.com.

[Giovanna Tramontin]

Dear Professor Fang, 

thank you for your response. 

One more question, I want to select the cwcata and the faces corresponding to z = 10, and map them in terms of x and y. For instance, the cwdata is in terms of the faces, the faces in terms of nodes and the nodes in terms of x,y, z. Is there a function to map the cwdata in terms of the nodes so I can plot the cwdata in xy for z = 10? 

Best wishes,

Giovanna 

[Qianqian Fang]

you can use the meshcentroid() function in iso2mesh to compute the centroid of each surface triangle, and associate the dref value to that centroid (but you will have to tessellate these centroid to form a new mesh).

alternatively, you can split the dref value on each triangle and attribute 1/3*area*elem_value to each vertex of the triangle, sum them for each surface node, and then divide the nodevolume() of each node, which is the sum of the 1/3*area of all surrounding triangles attached to a node. pseudo code can be found below:

% surftri is the surface triangle node list


nodal_dref=zeros(1:length(cfg.node))

tri_area=elemvolume(cfg.node, surftri)

nodal_area=nodevolume(cfg.node, surftri, tri_area)

for i=1:length(surftri)

  for j=1:3

    nodal_dref(surftri(i,j)) = nodal_dref(surftri(i,j)) + tri_area(i)*dref(i)

  end

end

nodal_dref=nodal_dref*(1/3);

nodal_dref(nodal_area > 0)=nodal_dref(nodal_area > 0)./nodal_area(nodal_area > 0)

you can do this more efficiently if you use vector operations (i.e. use accumarray() function)

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mmc-users/0fec4778-cc45-4fa1-94a3-018b3543d1aan%40googlegroups.com.