add demo script for mirror bc

mcxlab/examples/demo_focus_mirror_bc.m

@@ -0,0 +1,33 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% MCXLAB - Monte Carlo eXtreme for MATLAB/Octave by Qianqina Fang
+%
+% In this demo script, we verify the following two bugs are fixed
+%
+%  https://github.com/fangq/mcx/issues/103  (to handle focal point)
+%  https://github.com/fangq/mcx/issues/104  (to handle mirror bc)
+%
+% This file is part of Monte Carlo eXtreme (MCX) URL:http://mcx.sf.net
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+clear cfg
+cfg.nphoton=1e7;
+cfg.vol=uint8(ones(60,60,60));
+cfg.srctype='pattern';
+cfg.srcpos=[-40,-40, 0];
+cfg.srcparam1=[80 0 0 100];
+cfg.srcparam2=[0 80 0 100];
+cfg.srcdir=[0 0 1 30];
+cfg.issrcfrom0=1;
+cfg.srcpattern=zeros(100,100);
+cfg.srcpattern(51:end,51:end)=1;
+cfg.gpuid=1;
+cfg.autopilot=1;
+cfg.prop=[0 0 1 1;0.005 0.1 0 1];
+cfg.tstart=0;
+cfg.tend=5e-9;
+cfg.tstep=5e-09;
+cfg.bc='mm____';
+
+flux=mcxlab(cfg);
+
+mcxplotvol(log10(double(flux.data)))

src/mcx_core.cu

@@ -991,6 +991,10 @@ __device__ inline int launchnewphoton(MCXpos *p,MCXdir *v,MCXtime *f,float3* rv,
 		      break;
 		}
 	    }
+
+	    if(p->w<=gcfg->minenergy)
+	        continue;
+
             /**
              * If beam focus is set, determine the incident angle
              */
@@ -1470,7 +1474,7 @@ kernel void mcx_main_loop(uint media[],OutputType field[],float genergy[],uint n
 	  if(isreflect){
 	      if(gcfg->mediaformat<100)
 	          updateproperty<islabel>(&prop,mediaid,t); ///< optical property across the interface
-	      if(((isreflect && (isdet & 0xF)==0) || ((isdet & 0xF) & 0x1)) && ((isdet & 0xF)==bcMirror || n1!=((gcfg->mediaformat<100)? (prop.n):(gproperty[(mediaid>0 && gcfg->mediaformat>=100)?1:mediaid].w)))){
+	      if(((isreflect && (isdet & 0xF)==0) || (isdet & 0x1)) && ((isdet & 0xF)==bcMirror || n1!=((gcfg->mediaformat<100)? (prop.n):(gproperty[(mediaid>0 && gcfg->mediaformat>=100)?1:mediaid].w)))){
 	          float Rtotal=1.f;
 	          float cphi,sphi,stheta,ctheta,tmp0,tmp1;