Correction for Z Gaussian beam

[Tigrane Cantat-Moltrecht]

Hi all, 

wev'e been struggling for weeks because the Zgaussian sources would only launch around the upper vertical direction.

Based on the rotate_vector function of mcx, we reimplemented the vector rotation inside the zgaussian code:

in mmc_core.cl, for GPU usage:

------------------------

#if defined(__NVCC__) || defined(MCX_SRC_ZGAUSSIAN)
#ifdef __NVCC__
    } else if (GPU_PARAM(gcfg, srctype) == MCX_SRC_ZGAUSSIAN) {
#endif
            float ang, stheta, ctheta, sphi, cphi;
            ang = TWO_PI * rand_uniform01(ran); //next arimuth angle
            sphi = MCX_MATHFUN(sin)(ang);
            cphi = MCX_MATHFUN(cos)(ang);
            ang = MCX_MATHFUN(sqrt)(-2.f * MCX_MATHFUN(log)((rand_uniform01(ran)))) * (1.f - 2.f * rand_uniform01(ran)) * gcfg->srcparam1.x;
            stheta = MCX_MATHFUN(sin)(ang);
            ctheta = MCX_MATHFUN(cos)(ang);

            float x_old, y_old, z_old, tmp0, tmp1;
            x_old = r->vec.x;
            y_old = r->vec.y;
            z_old = r->vec.z;

            if (z_old > -1.f + EPS && z_old < 1.f - EPS) {
                tmp0 = 1.f - z_old * z_old;
                tmp1 = stheta * MCX_MATHFUN(rsqrtf)(tmp0);
                r->vec.x = tmp1 * (x_old * z_old * cphi - y_old * sphi) + x_old * ctheta;
                r->vec.y = tmp1 * (y_old * z_old * cphi + x_old * sphi) + y_old * ctheta;
                r->vec.z = -tmp1 * tmp0 * cphi + z_old * ctheta;
            }
            else {
                r->vec.x = stheta * cphi;
                r->vec.y = stheta * sphi;
                r->vec.z = (z_old > 0.f) ? ctheta : -ctheta;
            }
            canfocus = 0;
#endif

-------------------

in mmc_raytrace.c for CPU usage:

-------------------

    } else if (cfg->srctype == stZGaussian) {
        float ang, stheta, ctheta, sphi, cphi;
        ang = TWO_PI * rand_uniform01(ran); //next azimuth angle
        sphi = sinf(ang);
        cphi = cosf(ang);
        ang = sqrtf(-2.f * log(rand_uniform01(ran))) * (1.f - 2.f * rand_uniform01(ran0)) * cfg->srcparam1.x;
        stheta = sinf(ang);
        ctheta = cosf(ang);

        float x_old, y_old, z_old, tmp0, tmp1;
        x_old = r->vec.x;
        y_old = r->vec.y;
        z_old = r->vec.z;

        if (z_old > -1.f + EPS && z_old < 1.f - EPS) {
            tmp0 = 1.f - z_old * z_old;
            tmp1 = stheta / sqrtf(tmp0);
            r->vec.x = tmp1 * (x_old * z_old * cphi - y_old * sphi) + x_old * ctheta;
            r->vec.y = tmp1 * (y_old * z_old * cphi + x_old * sphi) + y_old * ctheta;
            r->vec.z = -tmp1 * tmp0 * cphi + z_old * ctheta;
        }
        else {
            r->vec.x = stheta * cphi;
            r->vec.y = stheta * sphi;
            r->vec.z = (z_old > 0.f) ? ctheta : -ctheta;
        }
        canfocus = 0;

-------------------

In principle, this sould not be limited to Zgaussian and should also apply to disk source, etc. When possible I'll try to propose a more geenral update.

I also want to try to input an "arbitrary" source to simulate LEDs, with a theta angle drawn from a specsheet radiation pattern or something ...

Best,

Tigrane