SourceFunction.hpp

Go to the documentation of this file.

/**
 * \file SourceFunction.hpp
 * \example mofem/tutorials/scl-10/src/SourceFunction.hpp
 *
 * Source function for photon diffusion problem.
 */
 
#include <boost/math/quadrature/gauss_kronrod.hpp>
using namespace boost::math::quadrature;
 
#include <stdlib.h>
#include <cmath>
#include <MoFEM.hpp>
 
using namespace MoFEM;
 
 
namespace SourceFunction {
    /** 
     * @brief Pulse is infinitely short.
     *
     * \note It is good approximation of pulse in femtosecond scale. To make it
     * longer one can apply third integral over time.
     *
     * \note Note analysis in photon_diffusion is shifted in time by initial_time.
     *
     * @param x
     * @param y
     * @param z
     * @return double
     */
 
    //! [sourceFunction]
    double sourceFunctionEval(const double x, const double y, const double z, const double beam_radius, const double beam_centre_x, 
                                    const double beam_centre_y, const double slab_thickness, const double mu_a, const double mu_sp,
                                    const double flux_magnitude, double initial_time, const double v, const double D) {
        const double A = 4. * D * v * initial_time;
        const double T =
            (v / pow(M_PI * A, 3. / 2.)) * exp(-mu_a * v * initial_time);
 
        auto phi_pulse = [&](const double r_s, const double phi_s) {
            const double xs = r_s * cos(phi_s);
            const double ys = r_s * sin(phi_s);
            const double xp = x - xs - beam_centre_x;
            const double yp = y - ys - beam_centre_y;
            const double zp1 = z + slab_thickness / 2. - 1. / mu_sp;
            const double zp2 = z + slab_thickness / 2. + 1. / mu_sp;
            const double P1 = xp * xp + yp * yp + zp1 * zp1;
            const double P2 = xp * xp + yp * yp + zp2 * zp2;
            return r_s * (exp(-P1 / A) - exp(-P2 / A));
        };
 
        auto f = [&](const double r_s) {
            auto g = [&](const double phi_s) { return phi_pulse(r_s, phi_s); };
            return gauss_kronrod<double, 15>::integrate(
                g, 0, 2 * M_PI, 0, std::numeric_limits<float>::epsilon());
        };
 
        return T * flux_magnitude *
                gauss_kronrod<double, 15>::integrate(
                    f, 0, beam_radius, 0, std::numeric_limits<float>::epsilon());
    };
  //! [sourceFunction]
    };