fix cyclic bc, add mcxlab demo

Author: fangq

mcxlab/examples/demo_infinite_slab_cyclic_bc.m

@@ -0,0 +1,34 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% MCXLAB - Monte Carlo eXtreme for MATLAB/Octave by Qianqina Fang
+%
+% In this example, we show how to simulate an infinite slab using cyclic
+% boundary condition
+%
+% This file is part of Monte Carlo eXtreme (MCX) URL:http://mcx.space
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+clear cfg;
+cfg.nphoton=1e8;
+cfg.issrcfrom0=1;
+cfg.vol=uint8(ones(60,60,20));
+cfg.srcdir=[0 0 1];
+cfg.gpuid=1;
+cfg.autopilot=1;
+cfg.prop=[0 0 1 1;0.005 1 0.8 1.37];
+cfg.tstart=0;
+cfg.seed=99999;
+
+% a uniform planar source outside the volume
+cfg.srctype='planar';
+cfg.srcpos=[0 0 0];
+cfg.srcparam1=[60 0 0 0];
+cfg.srcparam2=[0 60 0 0];
+cfg.tend=5e-9;
+cfg.tstep=5e-9;
+cfg.bc='ccrccr';
+flux=mcxlab(cfg);
+
+fcw=flux.data*cfg.tstep;
+imagesc(log10(abs(squeeze(fcw(:,30,:)))))
+axis equal; colorbar
+title('a uniform planar source incident along an infinite slab');

src/mcx_core.cu

@@ -1280,15 +1280,17 @@ kernel void mcx_main_loop(uint media[],OutputType field[],float genergy[],uint n
 	  /** launch new photon when exceed time window or moving from non-zero voxel to zero voxel without reflection */
           if((mediaid==0 && (((isdet & 0xF)==0 && (!gcfg->doreflect || (gcfg->doreflect && n1==gproperty[0].w))) || (isdet==bcAbsorb || isdet==bcCylic) )) || f.t>gcfg->twin1){
 	      if(isdet==bcCylic){
-		 if(flipdir==0)  p.x=mcx_nextafterf(roundf(p.x+((idx1d==OUTSIDE_VOLUME_MIN) ? gcfg->maxidx.x: -gcfg->maxidx.x)),(v.x > 0.f)-(v.x < 0.f));
-		 if(flipdir==1)  p.y=mcx_nextafterf(roundf(p.y+((idx1d==OUTSIDE_VOLUME_MIN) ? gcfg->maxidx.y: -gcfg->maxidx.y)),(v.y > 0.f)-(v.y < 0.f));
-		 if(flipdir==2)  p.z=mcx_nextafterf(roundf(p.z+((idx1d==OUTSIDE_VOLUME_MIN) ? gcfg->maxidx.z: -gcfg->maxidx.z)),(v.z > 0.f)-(v.z < 0.f));
-                 idx1d=(int(floorf(p.z))*gcfg->dimlen.y+int(floorf(p.y))*gcfg->dimlen.x+int(floorf(p.x)));
-	         mediaid=media[idx1d];
-	         isdet=mediaid & DET_MASK;  /** upper 16bit is the mask of the covered detector */
-	         mediaid &= MED_MASK;       /** lower 16bit is the medium index */
-                 GPUDEBUG(("Cylic boundary condition, moving photon in dir %d at %d flag, new pos=[%f %f %f]\n",flipdir,isdet,p.x,p.y,p.z));
-	         continue;
+                 if(flipdir==0)  p.x=mcx_nextafterf(roundf(p.x+((idx1d==OUTSIDE_VOLUME_MIN) ? gcfg->maxidx.x: -gcfg->maxidx.x)),(v.x > 0.f)-(v.x < 0.f));
+                 if(flipdir==1)  p.y=mcx_nextafterf(roundf(p.y+((idx1d==OUTSIDE_VOLUME_MIN) ? gcfg->maxidx.y: -gcfg->maxidx.y)),(v.y > 0.f)-(v.y < 0.f));
+                 if(flipdir==2)  p.z=mcx_nextafterf(roundf(p.z+((idx1d==OUTSIDE_VOLUME_MIN) ? gcfg->maxidx.z: -gcfg->maxidx.z)),(v.z > 0.f)-(v.z < 0.f));
+		 if(!(p.x<0.f||p.y<0.f||p.z<0.f||p.x>=gcfg->maxidx.x||p.y>=gcfg->maxidx.y||p.z>=gcfg->maxidx.z)){
+                     idx1d=(int(floorf(p.z))*gcfg->dimlen.y+int(floorf(p.y))*gcfg->dimlen.x+int(floorf(p.x)));
+	             mediaid=media[idx1d];
+	             isdet=mediaid & DET_MASK;  /** upper 16bit is the mask of the covered detector */
+	             mediaid &= MED_MASK;       /** lower 16bit is the medium index */
+                     GPUDEBUG(("Cylic boundary condition, moving photon in dir %d at %d flag, new pos=[%f %f %f]\n",flipdir,isdet,p.x,p.y,p.z));
+	             continue;
+		 }
 	      }
               GPUDEBUG(("direct relaunch at idx=[%d] mediaid=[%d], ref=[%d] bcflag=%d timegate=%d\n",idx1d,mediaid,gcfg->doreflect,isdet,f.t>gcfg->twin1));
 	      if(launchnewphoton<mcxsource>(&p,&v,&f,&rv,&prop,&idx1d,field,&mediaid,&w0,(mediaidold & DET_MASK),ppath,