revert the mua change made yesterday for dmmc, thanks to Shijie

mmclab/example/demo_mmcl_b1_b1d.m

@@ -92,16 +92,16 @@
 dets=[xi(:) 30.2*ones(size(xi(:))) yi(:)];
 
 hold on
-[c h2]=contourf(xi,yi, log10(max(squeeze(Phicpu*cfg.tstep),1e-8)), clines, 'k-' );
-contour(xi,yi,log10(abs(Phigpu*cfg.tstep)),clines,'r:')
-contour(xi,yi,log10(abs(squeeze(Phidcpu(1:end-1,31,1:end-1))'*cfg.tstep)),clines,'b--')
-contour(xi,yi,log10(abs(squeeze(Phidgpu(1:end-1,31,1:end-1))'*cfg.tstep)),clines,'y--')
+[c h2]=contourf(xi,yi, log10(max(squeeze(Phicpu*cfg.tstep),1e-8)), clines, 'k-', 'linewidth', 2);
+contour(xi,yi,log10(abs(Phigpu*cfg.tstep)),clines,'m:','linewidth',2)
+contour(xi,yi,log10(abs(squeeze(Phidcpu(1:end-1,31,1:end-1))'*cfg.tstep)),clines,'b--','linewidth',2)
+contour(xi,yi,log10(abs(squeeze(Phidgpu(1:end-1,31,1:end-1))'*cfg.tstep)),clines,'g--','linewidth',2)
 
 axis equal  
 set(gca,'xlim',[0.5 59.5])
 set(gca,'fontsize',20)
 xlabel('x (mm)')
 ylabel('z (mm)')
-legend('MMC','MMCL','D-MMC','D-MMCL')
+legend('B1-MMC','B1-MMCL','B1D-MMC','B1D-MMCL')
 legend boxoff;
 box on;

mmclab/example/demo_mmcl_b2_b2d.m

@@ -35,7 +35,7 @@
 %%-----------------------------------------------------------------
 clear cfg
 
-cfg.nphoton=3e6;
+cfg.nphoton=1e6;
 cfg.seed=1648335518;
 cfg.srcpos=[30,30.1,0];
 cfg.srcdir=[0 0 1];

src/mmc_cl_host.c

@@ -475,29 +475,16 @@ is more than what your have specified (%d), please use the -H option to specify
      fieldlen=(fieldlen>>cfg->nbuffer);
      field=realloc(field,sizeof(field[0])*fieldlen);
      if(cfg->exportfield){
-         if(cfg->method==rtBLBadouelGrid){
+         if(cfg->basisorder==0 || cfg->method==rtBLBadouelGrid){
              for(i=0;i<fieldlen;i++)
 	         cfg->exportfield[i]+=field[i];
 	 }else{
-             float ww;
-	     if(cfg->basisorder==0){
-		 for(i=0;i<cfg->maxgate;i++)
-		   for(j=0;j<mesh->ne;j++){
-		     ww=field[i*mesh->ne+j];
-		     if(cfg->outputtype!=otEnergy && cfg->outputtype!=otWP)
-                         ww*=(mesh->med[mesh->type[j]].mua<=0.f) ? 0.f : (1.f/mesh->med[mesh->type[j]].mua);
-                     cfg->exportfield[i*mesh->nn+j]+=ww;
-		   }
-	     }else{
-		 for(i=0;i<cfg->maxgate;i++)
-		   for(j=0;j<mesh->ne;j++){
-		     ww=field[i*mesh->ne+j]*0.25f;
-		     if(cfg->outputtype!=otEnergy && cfg->outputtype!=otWP)
-                         ww*=(mesh->med[mesh->type[j]].mua<=0.f) ? 0.f : (1.f/mesh->med[mesh->type[j]].mua);
-		     for(int k=0;k<mesh->elemlen;k++)
-			 cfg->exportfield[i*mesh->nn+mesh->elem[j*mesh->elemlen+k]-1]+=ww;
-		   }
-	     }
+             for(i=0;i<cfg->maxgate;i++)
+	       for(j=0;j<mesh->ne;j++){
+		 float ww=field[i*mesh->ne+j]*0.25f;
+	         for(int k=0;k<mesh->elemlen;k++)
+	             cfg->exportfield[i*mesh->nn+mesh->elem[j*mesh->elemlen+k]-1]+=ww;
+	       }
 	 }
      }
 

src/mmc_core.cl

@@ -506,6 +506,10 @@ float branchless_badouel_raytet(ray *r, __constant MCXParam *gcfg,__constant int
                     tshift=MIN( ((int)((r->photontimer-gcfg->tstart)*gcfg->Rtstep)), gcfg->maxgate-1 )*gcfg->framelen;
             {
 #ifndef MCX_SKIP_VOLUME
+	       if(prop.mua>0.f){
+	           if(gcfg->outputtype!=otEnergy && gcfg->outputtype!=otWP)
+                      ww/=prop.mua;
+	       }
   #ifndef USE_DMMC
                uint newidx=eid+tshift;
 	       r->oldidx=(r->oldidx==ID_UNDEFINED)? newidx: r->oldidx;

src/tettracing.c

@@ -1070,7 +1070,7 @@ float branchless_badouel_raytet(ray *r, raytracer *tracer, mcconfig *cfg, visito
 				if(newidx!=r->oldidx){
 				    tracer->mesh->weight[r->oldidx]+=r->oldweight;
 				    r->oldidx=newidx;
-				    r->oldweight=0.0;
+				    r->oldweight=w0*totalloss;
 				}else
 				    r->oldweight+=w0*totalloss;
 				w0*=segloss;