diff --git a/mcxcloud/frontend/index.html b/mcxcloud/frontend/index.html
index 30d1b2d7..54965e50 100644
--- a/mcxcloud/frontend/index.html
+++ b/mcxcloud/frontend/index.html
@@ -793,7 +793,9 @@
Backend
"zgaussian",
"line",
"slit",
- "pencilarray"
+ "pencilarray",
+ "hyperboloid",
+ "ring"
]
},
"Pos": {
diff --git a/mcxlab/examples/demo_mcxlab_srctype.m b/mcxlab/examples/demo_mcxlab_srctype.m
index 3ac4edcf..ab07e331 100644
--- a/mcxlab/examples/demo_mcxlab_srctype.m
+++ b/mcxlab/examples/demo_mcxlab_srctype.m
@@ -370,8 +370,62 @@
axis equal; colorbar
title('a gaussian beam source');
+
%% test group 5
+clear cfg;
+figure;
+cfg.nphoton=1e6;
+cfg.vol=uint8(ones(60,60,60));
+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;
+cfg.tend=5e-11;
+cfg.tstep=5e-11;
+
+% a ring beam
+cfg.srctype='ring';
+cfg.srcpos=[30 30 -10];
+cfg.srcdir=[0 0 1];
+cfg.srcparam1=[15 10 0 0];
+cfg.srcparam2=[0 0 0 0];
+flux=mcxlab(cfg);
+fcw=flux.data*cfg.tstep;
+subplot(221);
+imagesc(log10(abs(squeeze(fcw(:,:,1)))))
+axis equal; colorbar
+title('a ring beam');
+
+% a diverging ring sector beam
+cfg.srctype='ring';
+cfg.srcpos=[30 30 -10];
+cfg.srcdir=[0 0 1 -20];
+cfg.srcparam1=[15 10 0 2*pi/3];
+cfg.srcparam2=[0 0 0 0];
+flux=mcxlab(cfg);
+fcw=flux.data*cfg.tstep;
+subplot(222);
+imagesc(log10(abs(squeeze(fcw(:,:,1)))))
+axis equal; colorbar
+title('a ring-sector beam');
+
+% a ring beam
+cfg.srctype='pencilarray';
+cfg.srcpos=[10 10 0];
+cfg.srcdir=[0 0 1];
+cfg.srcparam1=[40 0 0 4];
+cfg.srcparam2=[0 40 0 5];
+flux=mcxlab(cfg);
+fcw=flux.data*cfg.tstep;
+subplot(223);
+imagesc(log10(abs(squeeze(fcw(:,:,1)))))
+axis equal; colorbar
+title('a pencil beam array)');
+
+%% test group 6
+
%debug flag to retrieve/test build-in RNG
cfg.vol=uint8(ones(100,100,100));
cfg.debuglevel='R';
diff --git a/mcxstudio/mcxsource.lfm b/mcxstudio/mcxsource.lfm
index 8a815dff..402e1b11 100644
--- a/mcxstudio/mcxsource.lfm
+++ b/mcxstudio/mcxsource.lfm
@@ -45,6 +45,8 @@ object fmSource: TfmSource
'line'
'slit'
'pencilarray'
+ 'hyperboloid'
+ 'ring'
)
OnEditingDone = edSourceEditingDone
ParentFont = False
diff --git a/pymcx b/pymcx
index e957869b..c5c841d0 160000
--- a/pymcx
+++ b/pymcx
@@ -1 +1 @@
-Subproject commit e957869b5761d76e12a8e7e665171946828c2469
+Subproject commit c5c841d02e1252c83f88c73337d6c870bb8b668c
diff --git a/schema/mcxinput.json b/schema/mcxinput.json
index aa34f0fd..705df7f1 100644
--- a/schema/mcxinput.json
+++ b/schema/mcxinput.json
@@ -216,7 +216,9 @@
"zgaussian",
"line",
"slit",
- "pencilarray"
+ "pencilarray",
+ "hyperboloid",
+ "ring"
]
},
"Pos": {
diff --git a/src/mcx_const.h b/src/mcx_const.h
index 0cd6407c..7fbdb5fb 100644
--- a/src/mcx_const.h
+++ b/src/mcx_const.h
@@ -93,6 +93,7 @@
#define MCX_SRC_PENCILARRAY 14 /**< a rectangular array of pencil beams */
#define MCX_SRC_PATTERN3D 15 /**< a 3D pattern source, starting from srcpos, srcparam1.{x,y,z} define the x/y/z dimensions */
#define MCX_SRC_HYPERBOLOID_GAUSSIAN 16 /**< Gaussian-beam with spot focus, scrparam1.{x,y,z} define beam waist, distance from source to focus, rayleigh range */
+#define MCX_SRC_RING 17 /**< ring/ring-sector source, scrparam1.{x,y} defines the outer/inner radius, srcparam1.{z,w} defines start/end angle*/
#define SAVE_DETID(a) ((a) & 0x1) /**< mask to save detector ID*/
#define SAVE_NSCAT(a) ((a)>>1 & 0x1) /**< output partial scattering counts */
diff --git a/src/mcx_core.cu b/src/mcx_core.cu
index 3f4a023b..6fb78772 100644
--- a/src/mcx_core.cu
+++ b/src/mcx_core.cu
@@ -328,7 +328,7 @@ __device__ inline void updatestokes(Stokes* s, float theta, float phi, float3* u
temp = (u2->z > -1.f && u2->z < 1.f) ? rsqrtf((1.f - costheta * costheta) * (1.f - u2->z * u2->z)) : 0.f;
cosi = (temp == 0.f) ? 0.f : (((phi > ONE_PI && phi < TWO_PI) ? 1.f : -1.f) * (u2->z * costheta - u->z) * temp);
- cosi = fmax(-1.f, fmin(cosi, 1.f));
+ cosi = fmaxf(-1.f, fminf(cosi, 1.f));
sini = sqrtf(1.f - cosi * cosi);
cos22 = 2.f * cosi * cosi - 1.f;
@@ -1240,15 +1240,20 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime*
}
case (MCX_SRC_DISK):
+ case (MCX_SRC_RING):
case (MCX_SRC_GAUSSIAN): { // uniform disk distribution or collimated Gaussian-beam
// Uniform disk point picking
// http://mathworld.wolfram.com/DiskPointPicking.html
- float sphi, cphi;
- float phi = TWO_PI * rand_uniform01(t);
+ float phi, sphi, cphi, r;
+
+ if (gcfg->srctype == MCX_SRC_RING && (gcfg->srcparam1.z > 0.f || gcfg->srcparam1.w > 0.f)) {
+ phi = fabsf(gcfg->srcparam1.z - gcfg->srcparam1.w) * rand_uniform01(t) + fminf(gcfg->srcparam1.z, gcfg->srcparam1.w);
+ } else {
+ phi = TWO_PI * rand_uniform01(t);
+ }
sincosf(phi, &sphi, &cphi);
- float r;
- if (gcfg->srctype == MCX_SRC_DISK) {
+ if (gcfg->srctype == MCX_SRC_DISK || gcfg->srctype == MCX_SRC_RING) {
r = sqrtf(rand_uniform01(t) * fabsf(gcfg->srcparam1.x * gcfg->srcparam1.x - gcfg->srcparam1.y * gcfg->srcparam1.y) + gcfg->srcparam1.y * gcfg->srcparam1.y);
} else if (fabsf(gcfg->c0.w) < 1e-5f || fabsf(gcfg->srcparam1.y) < 1e-5f) {
r = sqrtf(0.5f * rand_next_scatlen(t)) * gcfg->srcparam1.x;
@@ -1724,7 +1729,7 @@ __global__ void mcx_main_loop(uint media[], OutputType field[], float genergy[],
// in early CUDA, when ran=1, CUDA gives 1.000002 for tmp0 which produces nan later
// detected by Ocelot,thanks to Greg Diamos,see http://bit.ly/cR2NMP
- tmp0 = fmax(-1.f, fmin(1.f, tmp0));
+ tmp0 = fmaxf(-1.f, fminf(1.f, tmp0));
theta = acosf(tmp0);
stheta = sinf(theta);
@@ -1854,7 +1859,7 @@ __global__ void mcx_main_loop(uint media[], OutputType field[], float genergy[],
GPUDEBUG(("p=[%f %f %f] -> <%f %f %f>*%f -> hit=[%d %d %d] flip=%d\n", p.x, p.y, p.z, v.x, v.y, v.z, len, flipdir[0], flipdir[1], flipdir[2], flipdir[3]));
/** if the consumed unitless scat length is less than what's left in f.pscat, keep moving; otherwise, stop in this voxel */
- slen = fmin(slen, f.pscat);
+ slen = fminf(slen, f.pscat);
/** final length that the photon moves - either the length to move to the next voxel, or the remaining scattering length */
len = ((prop.mus == 0.f) ? len : (slen / prop.mus * (v.nscat + 1.f > gcfg->gscatter ? (1.f - prop.g) : 1.f)));
diff --git a/src/mcx_utils.c b/src/mcx_utils.c
index 05ad4c40..8c3aec07 100644
--- a/src/mcx_utils.c
+++ b/src/mcx_utils.c
@@ -196,7 +196,7 @@ const char boundarydetflag[] = {'0', '1', '\0'};
const char* srctypeid[] = {"pencil", "isotropic", "cone", "gaussian", "planar",
"pattern", "fourier", "arcsine", "disk", "fourierx", "fourierx2d", "zgaussian",
- "line", "slit", "pencilarray", "pattern3d", "hyperboloid", ""
+ "line", "slit", "pencilarray", "pattern3d", "hyperboloid", "ring", ""
};
diff --git a/src/mcx_utils.h b/src/mcx_utils.h
index 018424ff..ddcf44f4 100644
--- a/src/mcx_utils.h
+++ b/src/mcx_utils.h
@@ -211,7 +211,8 @@ typedef struct MCXConfig {
char internalsrc; /**<1 all photons launch positions are inside non-zero voxels, 0 let mcx search entry point*/
int zipid; /** 0)
+ [ncyl,fcyl]=meshacylinder(srcpos,srcpos+cfg(i).srcdir(1:3)*1e-5,cfg(i).srcparam1(2)*voxelsize,0,0);
+ hsrcarea=plotmesh(ncyl,fcyl{end-1},'facecolor','k','linestyle','none');
+ end
elseif(strcmp(cfg(i).srctype,'planar') || strcmp(cfg(i).srctype,'pattern') || strcmp(cfg(i).srctype,'fourier') || ...
strcmp(cfg(i).srctype,'fourierx') || strcmp(cfg(i).srctype,'fourierx2d') || strcmp(cfg(i).srctype,'pencilarray'))
if(~isfield(cfg(i),'srcparam1') || ~isfield(cfg(i),'srcparam2'))