Incorrect reflection (nothing reflects for small Fresnel coefficient)

[devhello145]

When Fresnel coefficient is small (< 0.02), photons are not reflected, even with large number nphoton (1e+7). Maybe the problem is in random number generator.

---

In this example of code, I create a 2layer mesh (1 - layer is transparent with n1 refractive index, 2-layer is absorb everything with n2 refractive index). The source is pencil, place on the botton of layer1 and incident angle (alpha) is const. So we know position of incident on boundary of layer 1 and layer 2.

The detector is placed in the position of reflected photons (in bottom of layer 1).

Also I add Fresnel calculation coefficient (r) for MMC (take from source code and slightly modify for octave) and MCML (see link), but they are almost same.

So we expect to detect (r * nphoton), when nphoton is enough large. But It does not (detect 0 photons)!!!WHY???.

---

See MAIN PARAMS section for setup:

• when n1 = 1.3 (we detect nothing!!!! for CPU and GPU mode) WHY???
• when n1 = 1.2 (we detect, but ratio of detected to nphoton is smaller then r). Maybe this another bug.

%% add path to isotomesh and mmclab

%% prepare mesh
h = 1;
[node, face, c0] = latticegrid([-20 20], [-20 20], [0, h, h + 10]);
% surf2mesh: v,f,p0,p1,keepratio,maxvol,regions,holes,forcebox
[node, elem] = surf2mesh(node, face, [], [], 1, 1, c0);
% plotmesh(node, elem);

%% MAIN PARAMS
srcpos = [0. 0. 0.];
alpha = 55; % incident angle
srcdir = [sind(alpha) 0 sqrt(1 - sind(alpha)^2)];
detpos = [2 * tand(alpha) * h, 0, 0, 5]; % set reflected array

n1 = 1.3;
n2 = 1.5;
prop = [0, 0, 1, 1; 0, 0, 1, n1; 1e+4, 1e-6, 1, n2];

nphoton = 1e+7;
gpuid = -1; % -1 = CPU, 1 = GPU

%% Functions for calculation of Fresnel reflectance
function r = reflectance_mmc(alpha, n1, n2)
    %% see "reflectray" function in file "mmc_raytrace.c"
    %% this code is slightly modified for octave

    Icos = cosd(alpha);

    tmp0 = n1 * n1;
    tmp1 = n2 * n2;
    tmp2 = 1 - tmp0 / tmp1 * (1 .- Icos * Icos); %*1-[n1/n2*sin(si)]^2 = cos(ti)^2*/

    % if (tmp2 > 0.f && !(*eid <= 0 && cfg->isreflect == bcMirror)) { /*if no total internal reflection*/
    if tmp2 > 0
        Re = tmp0 * Icos * Icos + tmp1 * tmp2; %/*transmission angle*/
        tmp2 = sqrt(tmp2); %/*to save one sqrt*/
        Im = 2 * n1 * n2 * Icos * tmp2;
        Rtotal = (Re - Im) / (Re + Im); %/*Rp*/
        Re = tmp1 * Icos * Icos + tmp0 * tmp2 * tmp2;
        Rtotal = (Rtotal + (Re - Im) / (Re + Im)) * 0.5; %/*(Rp+Rs)/2*/
    else
        Rtotal = 1;
    end
    r = Rtotal;
end

function r = reflectance_mcml(alpha, n1, n2)
    %% Take from MCML code, see "RFresnel" function in file "mcmlgo.c"
    ca1 = cosd(alpha);

    sa1 = sqrt(1 - ca1 * ca1);
    sa2 = n1 * sa1 / n2;
    if (sa2 >= 1.0)
        %/* double check for total internal reflection. */
        %*ca2_Ptr = 0.0;
        r = 1.0;
    else
        %   double cap, cam;    /* cosines of the sum ap or */
        %                   /* difference am of the two */
        %                   /* angles. ap = a1+a2 */
        %                   /* am = a1 - a2. */
        %   double sap, sam;    /* sines. */

        %   *ca2_Ptr = ca2 = sqrt(1-sa2*sa2);
        ca2 = sqrt(1 - sa2 * sa2);

        cap = ca1 * ca2 - sa1 * sa2; %/* c+ = cc - ss. */
        cam = ca1 * ca2 + sa1 * sa2; %/* c- = cc + ss. */
        sap = sa1 * ca2 + ca1 * sa2; %/* s+ = sc + cs. */
        sam = sa1 * ca2 - ca1 * sa2; %/* s- = sc - cs. */
        r = 0.5 * sam * sam * (cam * cam + cap * cap) / (sap * sap * cam * cam);
        % /* rearranged for speed. */
    end
end

%% CONFIG
cfg = struct();
cfg.seed = 1; %1648335518;

cfg.srctype = 'pencil';
cfg.srcpos = srcpos;
cfg.srcdir = srcdir;

%% set time of flight  (set to maximum)
cfg.tstart = 0;
cfg.tend = 1e-6;
cfg.tstep = 1e-6;

cfg.isreflect = 1;
cfg.issaveexit = 1;
cfg.issaveseed = 1; % overrides in mmclab by number of outputs

cfg.outputtype = 'fluence';

%% pernament
cfg.unitinmm = 1.;

%% tuned
cfg.minenergy = 1e-9;
cfg.maxdetphoton = 1e+7; % 1e+8 - broken
cfg.method = 'elem';
cfg.gpuid = -1; % -1 - CPU mode, 1 - GPU mode

%% additional
cfg.debuglevel = 'TP'; % M - movement info;

cfg.prop = prop;
cfg.detpos = detpos;

%% pernament (mesh)
cfg.node = node;
cfg.elem = elem(:, 1:4);
cfg.elemprop = elem(:, 5);

%% tuned
cfg.nphoton = nphoton;
cfg.gpuid = gpuid;

%% RUN
[fluence, detp, ncfg, seeds] = mmclab(cfg);

%% POST PROCESSING
r_mmc = reflectance_mmc(alpha, n1, n2);
r_mcml = reflectance_mcml(alpha, n1, n2);
fprintf('Fresnel Reflectance coefficients:\n\tMMC=  %.4f\n\tMCML= %.4f\n', r_mmc, r_mcml);
fprintf('Should be detected approximately: %.0f photons\n', r_mmc * cfg.nphoton);

if isempty(detp.data)
    fprintf('\nDetected 0 photons!\n\n');
    return
end
N = size(detp.ppath, 1);
fprintf('Detected: %i\n', N);

mask = detp.ppath(:, 2) == 0;
fprintf('Good: %i\n', sum(mask));
fprintf('Ratio: %f\n', sum(mask) / ncfg.nphoton);

---

Output for n1=1.3:

Cpu (gpuid=-1)//Nothing detects!!!:

```txt generating tetrahedral mesh from closed surfaces ... creating volumetric mesh from a surface mesh ... volume mesh generation is complete Launching MMCLAB - Mesh-based Monte Carlo for MATLAB & GNU Octave ... Running simulations for configuration #1 ... mmc.seed=1; mmc.srctype='pencil'; mmc.srcpos=[0 0 0]; mmc.srcdir=[0.819152 0 0.573576 0]; mmc.tstart=0; mmc.tend=1e-06; mmc.tstep=1e-06; mmc.isreflect=1; mmc.issaveexit=1; mmc.issaveseed=1; mmc.outputtype='fluence'; mmc.unitinmm=1; mmc.minenergy=1e-09; mmc.maxdetphoton=1e+07; mmc.method='elem'; mmc.gpuid=-1; mmc.debuglevel='TP'; mmc.prop=2; mmc.detnum=1; mmc.nn=7062; mmc.elem=[36761,4]; mmc.ne=36761; mmc.nphoton=1e+07; mmc.facenb=[36761,4]; mmc.evol=36761; mmc.e0=1874; done 3 simulating ... ############################################################################### # Mesh-based Monte Carlo (MMC) - OpenCL # # Copyright (c) 2010-2020 Qianqian Fang # # http://mcx.space/#mmc # # # #Computational Optics & Translational Imaging (COTI) Lab [http://fanglab.org]# # Department of Bioengineering, Northeastern University, Boston, MA, USA # # # # Research funded by NIH/NIGMS grant R01-GM114365 # ############################################################################### $Rev:: $v2021.2$Date:: $ by $Author:: $ ############################################################################### seed=1 simulating ... Progress: [==============================================================] 100% done 98264 speed ... 101.77 photon/ms, 11588222976 ray-tetrahedron tests (0 overhead, 117929.48 test/ms) detected 0 photons source 1 total simulated energy: 10000000.000000 absorbed: 100.00000% normalizor=1.38933e-07 saving fluence ...saving detected photons ... done 98266 done 98269 Fresnel Reflectance coefficients: MMC= 0.0154 MCML= 0.0154 Should be detected approximately: 153516 photons Detected 0 photons! ```

GPU (gpuid=1)//nothing detects!!!

```txt generating tetrahedral mesh from closed surfaces ... creating volumetric mesh from a surface mesh ... volume mesh generation is complete Launching MMCLAB - Mesh-based Monte Carlo for MATLAB & GNU Octave ... Running simulations for configuration #1 ... mmc.seed=1; mmc.srctype='pencil'; mmc.srcpos=[0 0 0]; mmc.srcdir=[0.819152 0 0.573576 0]; mmc.tstart=0; mmc.tend=1e-06; mmc.tstep=1e-06; mmc.isreflect=1; mmc.issaveexit=1; mmc.issaveseed=1; mmc.outputtype='fluence'; mmc.unitinmm=1; mmc.minenergy=1e-09; mmc.maxdetphoton=1e+07; mmc.method='elem'; mmc.gpuid=1; mmc.debuglevel='TP'; mmc.prop=2; mmc.detnum=1; mmc.nn=7062; mmc.elem=[36761,4]; mmc.ne=36761; mmc.nphoton=1e+07; mmc.facenb=[36761,4]; mmc.evol=36761; mmc.e0=1874; done 4 simulating ... ############################################################################### # Mesh-based Monte Carlo (MMC) - OpenCL # # Copyright (c) 2010-2020 Qianqian Fang # # http://mcx.space/#mmc # # # #Computational Optics & Translational Imaging (COTI) Lab [http://fanglab.org]# # Department of Bioengineering, Northeastern University, Boston, MA, USA # # # # Research funded by NIH/NIGMS grant R01-GM114365 # ############################################################################### $Rev:: $v2021.2$Date:: $ by $Author:: $ ############################################################################### - variant name: [MMC-OpenCL] compiled with OpenCL version [1] - compiled with: [RNG] xorshift128+ RNG [Seed Length] 4 initializing streams ... init complete : 0 ms Building kernel with option: -cl-mad-enable -DMCX_USE_NATIVE -DMCX_SIMPLIFY_BRANCH -DMCX_VECTOR_INDEX -DMCX_SRC_PENCIL -DUSE_ATOMIC -DMCX_SAVE_DETECTORS -DMCX_SAVE_SEED -DMCX_DO_REFLECTION -DUSE_BLBADOUEL build program complete : 3 ms - [device 0(1): NVIDIA GeForce RTX 2060 SUPER] threadph=71 oddphotons=112256 np=10000000.0 nthread=139264 nblock=64 repetition=1 set kernel arguments complete : 3 ms 3 lauching mcx_main_loop for time window [0.0ns 1000.0ns] ... simulation run# 1 ... Progress: [==============================================================] 100% kernel complete: 8550 ms retrieving flux ... detected 0 photons, total: 0 transfer complete: 8750 ms normalizing raw data ... simulated 10000000 photons (10000000) with 1 devices (ray-tet 11582557184) MCX simulation speed: 1170.00 photon/ms total simulated energy: 10000000.00 absorbed: 100.00000% (loss due to initial specular reflection is excluded in the total) done 9267 Fresnel Reflectance coefficients: MMC= 0.0154 MCML= 0.0154 Should be detected approximately: 153516 photons Detected 0 photons! ```

Output for n1=1.2: (please modify code example: set n1)

GPU (gpuid=1) // Detects smaller and not all good photons (but it's maybe another bag)

```txt generating tetrahedral mesh from closed surfaces ... creating volumetric mesh from a surface mesh ... volume mesh generation is complete Launching MMCLAB - Mesh-based Monte Carlo for MATLAB & GNU Octave ... Running simulations for configuration #1 ... mmc.seed=1; mmc.srctype='pencil'; mmc.srcpos=[0 0 0]; mmc.srcdir=[0.819152 0 0.573576 0]; mmc.tstart=0; mmc.tend=1e-06; mmc.tstep=1e-06; mmc.isreflect=1; mmc.issaveexit=1; mmc.issaveseed=1; mmc.outputtype='fluence'; mmc.unitinmm=1; mmc.minenergy=1e-09; mmc.maxdetphoton=1e+07; mmc.method='elem'; mmc.gpuid=1; mmc.debuglevel='TP'; mmc.prop=2; mmc.detnum=1; mmc.nn=7062; mmc.elem=[36761,4]; mmc.ne=36761; mmc.nphoton=1e+07; mmc.facenb=[36761,4]; mmc.evol=36761; mmc.e0=1874; done 3 simulating ... ############################################################################### # Mesh-based Monte Carlo (MMC) - OpenCL # # Copyright (c) 2010-2020 Qianqian Fang # # http://mcx.space/#mmc # # # #Computational Optics & Translational Imaging (COTI) Lab [http://fanglab.org]# # Department of Bioengineering, Northeastern University, Boston, MA, USA # # # # Research funded by NIH/NIGMS grant R01-GM114365 # ############################################################################### $Rev:: $v2021.2$Date:: $ by $Author:: $ ############################################################################### - variant name: [MMC-OpenCL] compiled with OpenCL version [1] - compiled with: [RNG] xorshift128+ RNG [Seed Length] 4 initializing streams ... init complete : 0 ms Building kernel with option: -cl-mad-enable -DMCX_USE_NATIVE -DMCX_SIMPLIFY_BRANCH -DMCX_VECTOR_INDEX -DMCX_SRC_PENCIL -DUSE_ATOMIC -DMCX_SAVE_DETECTORS -DMCX_SAVE_SEED -DMCX_DO_REFLECTION -DUSE_BLBADOUEL build program complete : 3 ms - [device 0(1): NVIDIA GeForce RTX 2060 SUPER] threadph=71 oddphotons=112256 np=10000000.0 nthread=139264 nblock=64 repetition=1 set kernel arguments complete : 3 ms 3 lauching mcx_main_loop for time window [0.0ns 1000.0ns] ... simulation run# 1 ... Progress: [==============================================================] 100% kernel complete: 589 ms retrieving flux ... detected 254771 photons, total: 254771 transfer complete: 792 ms normalizing raw data ... simulated 10000000 photons (10000000) with 1 devices (ray-tet 510743744) MCX simulation speed: 17064.85 photon/ms total simulated energy: 10000000.00 absorbed: 97.76932% (loss due to initial specular reflection is excluded in the total) done 1290 Fresnel Reflectance coefficients: MMC= 0.0302 MCML= 0.0302 Should be detected approximately: 301522 photons Detected: 254771 Good: 223059 Ratio: 0.022306 ```

[fangq]

mmc treats the space outside of the mesh (i.e. the entire space outside of your mesh) as medium 0 (background) and use cfg.prop(1,:)=[0,0,1,1] for such medium.

when you calculate your reflection ratio, did you consider the secondary reflection between your layer1 and the background medium?

[devhello145]

As I understand, there are only one reflection layer1 -> layer2. I took care for reflection of "layer2-> outside". It is impossible due to high mu_a and low mu_s of layer2. So all photons will die in layer2.

Is it possible that photons reflects from layer2, go to layer 1 and then reflects from outside?

---

when you calculate your reflection ratio, did you consider the secondary reflection between your layer1 and the background medium?

But I also experimented without secondary reflection by setting

n1=1.3;
prop = [0, 0, 1, n1; 0, 0, 1, n1; 1e+4, 1e-6, 1, n2];
gpuid = 1;

and detected nothing (0 photons on PD). Therefore the issue is still exist.

[fangq]

again, sorry for the long delay in getting back to you.

I ran your simulation and plotted a trajectory of a photon, see below

you can see that the photon is bounced at nearly all interfaces due to total internal reflection - with the exception that the y/z plane interface of medium type 1 (layer 1).

you can set your detpos to have a large radius (such as detpos = [2 * tand(alpha) * h, 0, 0, 50];) to capture all photons and plot the exiting position using

mcxpreview(cfg);
hold on
plotmesh(detp.p, 'bo')

you can see that all photons exit from the two strips along y=-20 and y=20.

from your domain setting, I believe mmc traces the photon properly and there should be no photon exiting from the z=0 plane.

I am closing this ticket. if you see any issue, feel free to reopen.