mccc_coefs.h

Go to the documentation of this file.

#ifndef MCCC_COEFS_H
#define MCCC_COEFS_H
 
#include <math.h>
#include "../../ascot5.h"
#include "../../consts.h"
#include "mccc.h"

#define mccc_coefs_cab(qa, qb, nb, clogab) (                            \
        nb * qa*qa * qb*qb * clogab / ( 4 * CONST_PI * CONST_E0*CONST_E0 ) )

#define mccc_coefs_Q(ma, qa, mb, qb, nb, vb, clogab, mu0) (           \
        -mccc_coefs_cab(qa, qb, nb, clogab) * mu0 / ( ma * mb * vb*vb ) )

#define mccc_coefs_dQ(ma, qa, mb, qb, nb, vb, clogab, dmu0) (           \
        -mccc_coefs_cab(qa, qb, nb, clogab) * dmu0 / ( ma * mb * vb*vb*vb ) )

#define mccc_coefs_F(ma, qa, mb, qb, nb, vb, clogab, mu0) (             \
        -( 1/mb + 1/ma ) * mccc_coefs_cab(qa, qb, nb, clogab) * mu0     \
        / ( ma * vb*vb ) )

#define mccc_coefs_Dpara(ma, qa, va, qb, nb, vb, clogab, mu0) (         \
        ( va > 0 ) ?                                                    \
        mccc_coefs_cab(qa, qb, nb, clogab) * mu0 / ( 2 * ma*ma * va ) : \
        mccc_coefs_cab(qa, qb, nb, clogab) * 4                          \
        / ( 6 * CONST_SQRTPI * ma*ma * vb ) )

#define mccc_coefs_dDpara(ma, qa, va, qb, nb, vb, clogab, mu0, dmu0) (  \
        mccc_coefs_cab(qa, qb, nb, clogab) * ( dmu0/vb - mu0/va )       \
        / ( 2 * ma*ma * va ) )

#define mccc_coefs_Dperp(ma, qa, va, qb, nb, vb, clogab, mu1) (         \
        ( va > 0 ) ?                                                    \
        mccc_coefs_cab(qa, qb, nb, clogab) * mu1 / ( 2 * ma*ma * va ) : \
        mccc_coefs_cab(qa, qb, nb, clogab) * 4                          \
        / ( 6 * CONST_SQRTPI * ma*ma * vb ) )

#define mccc_coefs_K(va, Dpara, dDpara, Q) (    \
        Q + dDpara + 2*Dpara / va )

 #define mccc_coefs_nu(va, Dperp) (             \
        2 * Dperp / ( va * va ) )

#define mccc_coefs_DX(xi, Dpara, Dperp, gyrofreq) (              \
        ( 0.5 * ( Dpara - Dperp ) * ( 1 - xi*xi ) + Dperp )      \
        / (gyrofreq*gyrofreq) )
 
GPU_DECLARE_TARGET_SIMD_UNIFORM(mdata)
static void mccc_coefs_mufun(real mufun[3], real x, mccc_data* mdata);
DECLARE_TARGET_END
 
GPU_DECLARE_TARGET_SIMD_UNIFORM(nspec, mb, qb, nb, Tb)
inline static void mccc_coefs_clog(
        real* clogab, real ma, real qa, real va, int nspec,
        const real* mb, const real* qb, const real* nb, const real* Tb);
DECLARE_TARGET_END

inline static void mccc_coefs_clog(
        real* clogab, real ma, real qa, real va, int nspec, const real* mb,
        const real* qb, const real* nb, const real* Tb) {
 
    /* Evaluate Debye length */
    real sum = 0;
    GPU_SEQUENTIAL_LOOP
    for(int i = 0; i < nspec; i++){
        sum += nb[i] * qb[i] * qb[i] / Tb[i];
    }
    real debyeLength = sqrt(CONST_E0/sum);
 
    /* Evaluate classical and quantum mechanical impact parameter. The one *
     * that is larger is used to evaluate Coulomb logarithm.               */
    GPU_SEQUENTIAL_LOOP
    for(int i=0; i < nspec; i++){
        real vbar = va * va + 2 * Tb[i] / mb[i];
        real mr   = ma * mb[i] / ( ma + mb[i] );
        real bcl  = fabs( qa * qb[i] / ( 4*CONST_PI*CONST_E0 * mr * vbar ) );
        real bqm  = fabs( CONST_HBAR / ( 2 * mr * sqrt( vbar ) ) );
 
        if(bcl > bqm){
            clogab[i] = log( debyeLength / bcl );
        }
        else{
            clogab[i] = log( debyeLength / bqm );
        }
    }
}

inline static void mccc_coefs_mufun(real mufun[3], real x, mccc_data* mdata) {
 
    if(!mdata->usetabulated && x!= 0) {
        real expm2x = exp(-x*x);
        real erfx   = erf(x);
 
        mufun[0] = ( erfx - 2 * x * expm2x / CONST_SQRTPI ) / (x*x);
        mufun[1] = erfx - 0.5 * mufun[0];
        mufun[2] = 4 * expm2x / CONST_SQRTPI - 2 * mufun[0] / x;
    }
    else if(mdata->usetabulated && x != 0) {
        // TODO implement me
    }
    else {
        mufun[0] = 0;
        mufun[1] = 0;
        mufun[2] = 4 / ( 3 * CONST_SQRTPI );
    }
 
}
#endif