normalisation.h

Go to the documentation of this file.

 
#ifndef NORMALISATION_H
#define NORMALISATION_H
 
#include "qpms_types.h"
#include "qpms_error.h"
#include <math.h>
#include <complex.h>
#include "indexing.h"
#include "optim.h"
 
 
static inline double qpms_normalisation_normfactor(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
    switch (QPMS_EXPECT(norm & QPMS_NORMALISATION_NORM_BITS, norm & QPMS_NORMALISATION_DEFAULT)) {
        case QPMS_NORMALISATION_NORM_POWER:
            return 1;
        case QPMS_NORMALISATION_NORM_SPHARM:
            return sqrt(l*(l+1));
        case QPMS_NORMALISATION_NORM_NONE: // TODO more precision
            return sqrt(l*(l+1) * 4*M_PI / (2*l+1)) *
                    exp(0.5*(lgamma(l+m+1) - lgamma(l-m+1)));
        default:
            QPMS_WTF;
    }
}
 
 
 
 
static inline complex double qpms_normalisation_factor_M_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
    complex double fac = qpms_normalisation_normfactor(norm, l, m);
    if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_M_MINUS)) fac *= -1;
    if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_M_I)) fac *= I;
    if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac;
    return fac;
}
 
 
 
static inline complex double qpms_normalisation_factor_M(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
    complex double fac = qpms_normalisation_factor_M_noCS(norm, l, m);
    return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac;
}
 
 
 
static inline complex double qpms_normalisation_factor_N_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
    complex double fac = qpms_normalisation_normfactor(norm, l, m);
    if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_N_MINUS)) fac *= -1;
    if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_N_I)) fac *= I;
    if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac;
    return fac;
}
 
 
 
static inline complex double qpms_normalisation_factor_N(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
    complex double fac = qpms_normalisation_factor_N_noCS(norm, l, m);
    return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac;
}
 
 
static inline complex double qpms_normalisation_factor_N_M(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
    return qpms_normalisation_factor_N_noCS(norm, l, m) 
        / qpms_normalisation_factor_M_noCS(norm, l, m);
}
 
 
 
static inline complex double qpms_normalisation_factor_L_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
    complex double fac = qpms_normalisation_normfactor(norm, l, m);
    if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_L_MINUS)) fac *= -1;
    if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_L_I)) fac *= I;
    if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac;
    return fac;
}
 
 
static inline complex double qpms_normalisation_factor_L(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
    complex double fac = qpms_normalisation_factor_L_noCS(norm, l, m);
    return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac;
}
 
static inline complex double qpms_normalisation_factor_uvswfi(const qpms_normalisation_t norm, qpms_uvswfi_t ui) {
    qpms_vswf_type_t t; qpms_m_t m; qpms_l_t l;
    qpms_uvswfi2tmn(ui, &t, &m, &l);
    switch(t) {
        case QPMS_VSWF_MAGNETIC:
            return qpms_normalisation_factor_M(norm, l, m);
        case QPMS_VSWF_ELECTRIC:
            return qpms_normalisation_factor_N(norm, l, m);
        case QPMS_VSWF_LONGITUDINAL:
            return qpms_normalisation_factor_L(norm, l, m);
        default:
            QPMS_WTF;
    }
}
 
 
 
static inline qpms_normalisation_t qpms_normalisation_dual(qpms_normalisation_t norm) {
    norm ^= QPMS_NORMALISATION_INVERSE;
    if (!(norm & QPMS_NORMALISATION_SPHARM_REAL))
        norm ^= QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE;
    return norm;
}
 
 
static inline complex double qpms_spharm_azimuthal_part(qpms_normalisation_t norm, qpms_m_t m, double phi) {
    switch(QPMS_EXPECT(norm, QPMS_NORMALISATION_DEFAULT) 
            & (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) {
        case 0:
            return cexp(I*m*phi);
        case QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE:
            return cexp(-I*m*phi);
        case QPMS_NORMALISATION_SPHARM_REAL:
            if (m > 0) return M_SQRT2 * cos(m*phi);
            else if (m < 0) return M_SQRT2 * sin(m*phi);
            else return 1.;
        default:
            QPMS_WTF;
    }
}
 
 
static inline complex double qpms_spharm_azimuthal_part_derivative_div_m(qpms_normalisation_t norm, qpms_m_t m, double phi) {
    if(m==0) return 0;
    switch(QPMS_EXPECT(norm, QPMS_NORMALISATION_DEFAULT) 
            & (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) {
        case 0:
            return I*cexp(I*m*phi);
        case QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE:
            return -I*cexp(-I*m*phi);
        case QPMS_NORMALISATION_SPHARM_REAL:
            if (m > 0) return -M_SQRT2 * sin(m*phi);
            else if (m < 0) return M_SQRT2 * cos(m*phi);
            else return -1;
        default:
            QPMS_WTF;
    }
}
 
#if 0 // legacy code moved from qpms_types.h. TODO cleanup
static inline int qpms_normalisation_t_normonly(qpms_normalisation_t norm) {
    return norm & (~QPMS_NORMALISATION_T_CSBIT);
}
 
/* Normalisation of the spherical waves is now scattered in at least three different files:
 * here, we have the norm in terms of radiated power of outgoing wave.
 * In file legendre.c, function qpms_pitau_get determines the norm used in the vswf.c
 * spherical vector wave norms. The "dual" waves in vswf.c use the ..._abssquare function below.
 * In file translations.c, the normalisations are again set by hand using the normfac and lognormfac
 * functions.
 */
#include <math.h>
#include <assert.h>
// relative to QPMS_NORMALISATION_KRISTENSSON_CS, i.e.
// P_l^m[normtype] = P_l^m[Kristensson]
static inline double qpms_normalisation_t_factor(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
    int csphase = qpms_normalisation_t_csphase(norm);
    norm = qpms_normalisation_t_normonly(norm);
    double factor;
    switch (norm) {
        case QPMS_NORMALISATION_KRISTENSSON:
            factor = 1.;
            break;
        case QPMS_NORMALISATION_TAYLOR:
            factor = sqrt(l*(l+1));
            break;
        case QPMS_NORMALISATION_NONE:
            factor = sqrt(l*(l+1) * 4 * M_PI / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1)));
            break;
#ifdef USE_XU_ANTINORMALISATION // broken probably in legendre.c
        case QPMS_NORMALISATION_XU:
            factor = sqrt(4 * M_PI) / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1));
            break;
#endif
        default:
            assert(0);
    }
    factor *= (m%2)?(-csphase):1;
    return factor;
}
 
 
// TODO move elsewhere
static inline double qpms_normalisation_t_factor_abssquare(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
    norm = qpms_normalisation_t_normonly(norm);
    switch (norm) {
        case QPMS_NORMALISATION_KRISTENSSON:
            return 1.;
            break;
        case QPMS_NORMALISATION_TAYLOR:
            return l*(l+1);
            break;
        case QPMS_NORMALISATION_NONE:
            return l*(l+1) * 4 * M_PI / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1));
            break;
#ifdef USE_XU_ANTINORMALISATION // broken probably in legendre.c
        case QPMS_NORMALISATION_XU:
            {
              double fac = sqrt(4 * M_PI) / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1));
              return fac * fac;
            }
            break;
#endif 
        default:
            assert(0);
            return NAN;
    }
}
#endif
 
#endif //NORMALISATION_H