Add gaussian broadening to SLIT source

- `param2.x`: width of gaussian broadening in direction perpendicular to both slit and v
- `param2.y`: width of gaussian broadening in direction of slit

src/mcx_core.cu

@@ -1462,6 +1462,29 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime*
                         rotatevector(v, 1.f, 0.f, sphi, cphi);
                     }
 
+                    if (MCX_SRC_SLIT && (launchsrc->param2.x > 0.f || launchsrc->param2.y > 0.f)) {
+                        float sphi, cphi;
+                        r = TWO_PI * rand_uniform01(t);
+                        sincosf(r, &sphi, &cphi);
+                        r = sqrtf(2.f * rand_next_scatlen(t));
+                        // gaussian broadening factor in direction perpendicular to both slit and v directions
+                        float gauss_perp = launchsrc->param2.x * r * cphi;
+                        // gaussian broadening factor in direction of slit
+                        float gauss_parallel = launchsrc->param2.y * r * sphi;
+                        float parallel_norm = rnorm3df(launchsrc->param1.x, launchsrc->param1.y, launchsrc->param1.z);
+                        float3 perp = float3(launchsrc->param1.y * v->z - launchsrc->param1.z * v->y,
+                                             launchsrc->param1.z * v->x - launchsrc->param1.x * v->z,
+                                             launchsrc->param1.x * v->y - launchsrc->param1.y * v->x);
+                        float perp_norm = rnorm3df(perp.x, perp.y, perp.z);
+                        v->x += gauss_perp * perp.x * perp_norm + gauss_parallel * launchsrc->param1.x * parallel_norm;
+                        v->y += gauss_perp * perp.y * perp_norm + gauss_parallel * launchsrc->param1.y * parallel_norm;
+                        v->z += gauss_perp * perp.z * perp_norm + gauss_parallel * launchsrc->param1.z * parallel_norm;
+                        r = rnorm3df(v->x, v->y, v->z);
+                        v->x *= r;
+                        v->y *= r;
+                        v->z *= r;
+                    }
+
                     *rv = float3(rv->x + (launchsrc->param1.x) * 0.5f,
                                  rv->y + (launchsrc->param1.y) * 0.5f,
                                  rv->z + (launchsrc->param1.z) * 0.5f);