157 lines
5.5 KiB
C
157 lines
5.5 KiB
C
/*! \file indexing.h
|
|
* \brief Various index conversion functions.
|
|
*
|
|
* Index conventions in QPMS
|
|
* -------------------------
|
|
*
|
|
* Depending on circumstances (e.g. whether scalar or vector spherical wavefunctions
|
|
* are being used), QPMS uses two mappings between the multipole degree/order index
|
|
* pair (\a l, \a m), due to historical reasons ocassionally also noted as (\a n, \a m),
|
|
* and a flat non-negative integer index.
|
|
*
|
|
* Both mappings map the (\a l, \a m) pair in lexicographic order (with
|
|
* the basic constraint \f$ |m| \le l \f$), the difference is that the "scalar"
|
|
* mapping starts from \a l = 0, whereas the "vector" mapping starts from \a l = 1,
|
|
* so that for the zeroth elements we have
|
|
* * (0, 0) ↔ 0, (1, -1) ↔ 1 in the "scalar" mapping,
|
|
* * (1, -1) ↔ 0 in the "vector" mapping.
|
|
*
|
|
* The vector SWF flat index is typically denoted by `y` and its type is \ref qpms_y_t,
|
|
* whereas * the scalar SWF flat index is typically denoted by `y_sc` and its
|
|
* type is \ref qpms_y_sc_t.
|
|
*
|
|
* Moreover, there is another mapping that encodes also the electric/magnetic/scalar
|
|
* WSWF type index, see \ref qpms_uvswfi_t.
|
|
*
|
|
*/
|
|
#ifndef QPMS_INDEXING_H
|
|
#define QPMS_INDEXING_H
|
|
#ifdef __cplusplus
|
|
extern "C" {
|
|
#endif
|
|
|
|
#include "qpms_types.h"
|
|
#include <math.h>
|
|
#include "optim.h"
|
|
|
|
/// Mapping from multipole degree, order to flat index (degree ≥ 1)
|
|
static inline qpms_y_t qpms_mn2y(
|
|
qpms_m_t m, ///< multipole order
|
|
qpms_l_t n ///< multipole degree (a.k.a. \a l); must be greater than 0
|
|
) {
|
|
return n * (n + 1) + m - 1;
|
|
}
|
|
|
|
/// Mapping from flat index to multipole degree (degree ≥ 1)
|
|
static inline qpms_lm_t qpms_y2n(qpms_y_t y) {
|
|
//return (sqrt(5+y)-2)/2; // the cast will truncate the fractional part, which is what we want
|
|
return sqrt(y+1);
|
|
}
|
|
|
|
/// Mapping from flat index, multipole degree (≥ 1) to multipole order
|
|
static inline qpms_m_t qpms_yn2m(qpms_y_t y, qpms_l_t n) {
|
|
return y-qpms_mn2y(0,n);
|
|
}
|
|
|
|
/// Mapping from flat index to multipole degree, order (degree ≥ 1)
|
|
static inline void qpms_y2mn_p(qpms_y_t y, qpms_m_t *m, qpms_l_t *n){
|
|
*m=qpms_yn2m(y,*n=qpms_y2n(y));
|
|
}
|
|
|
|
/// Number of spherical multipoles with `lmax` ≥ degree ≥ 1
|
|
static inline qpms_y_t qpms_lMax2nelem(qpms_l_t lmax){
|
|
return lmax * ((qpms_y_t)lmax + 2);
|
|
}
|
|
|
|
// Scalar versions: they have a place for the 0, 0 term in the beginning
|
|
|
|
/// Mapping from multipole degree, order to flat index (degree ≥ 0)
|
|
static inline qpms_y_sc_t qpms_mn2y_sc(
|
|
qpms_m_t m, ///< multipole order
|
|
qpms_l_t n ///< multipole degree (a.k.a. \a l); muste be greater or equal than 0
|
|
) {
|
|
return n * (n + 1) + m;
|
|
}
|
|
|
|
/// Mapping from flat index to multipole degree (degree ≥ 0)
|
|
static inline qpms_lm_t qpms_y2n_sc(qpms_y_sc_t y) {
|
|
//return (sqrt(5+y)-2)/2; // the cast will truncate the fractional part, which is what we want
|
|
return sqrt(y);
|
|
}
|
|
|
|
/// Mapping from flat index, multipole degree (≥ 0) to multipole order
|
|
static inline qpms_m_t qpms_yn2m_sc(qpms_y_sc_t y, qpms_l_t n) {
|
|
return y-qpms_mn2y_sc(0,n);
|
|
}
|
|
|
|
/// Mapping from flat index to multipole degree, order (degree ≥ 0)
|
|
static inline void qpms_y2mn_sc_p(qpms_y_sc_t y, qpms_m_t *m, qpms_l_t *n){
|
|
*m=qpms_yn2m_sc(y,*n=qpms_y2n_sc(y));
|
|
}
|
|
|
|
/// Number of spherical multipoles with `lmax` ≥ degree ≥ 0
|
|
static inline qpms_y_sc_t qpms_lMax2nelem_sc(qpms_l_t lmax){
|
|
return lmax * ((qpms_y_sc_t)lmax + 2) + 1;
|
|
}
|
|
|
|
// TODO maybe enable crashing / validity control by macro definitions...
|
|
|
|
/// Conversion from VSWF type, order and degree to universal index.
|
|
static inline qpms_uvswfi_t qpms_tmn2uvswfi(
|
|
qpms_vswf_type_t t, qpms_m_t m, qpms_l_t n) {
|
|
return t + 4 * qpms_mn2y_sc(m, n);
|
|
}
|
|
|
|
/// Constant denoting the l=0, m=0 longitudinal spherical vector wave.
|
|
static const qpms_uvswfi_t QPMS_UI_L00 = 0;
|
|
|
|
/// Conversion from universal VSWF index u to type, order and degree.
|
|
/** Returns a non-zero value if the u value is invalid. */
|
|
static inline qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u,
|
|
qpms_vswf_type_t *t, qpms_m_t *m, qpms_l_t *n) {
|
|
*t = (qpms_vswf_type_t)(u & 3);
|
|
qpms_y_sc_t y_sc = u / 4;
|
|
qpms_y2mn_sc_p(y_sc, m, n);
|
|
// Test validity
|
|
if (QPMS_UNLIKELY(*t == 3)) return QPMS_ERROR; // VSWF type code invalid, TODO WARN
|
|
if (QPMS_UNLIKELY(*t && !y_sc)) return QPMS_ERROR; // l == 0 for transversal wave, TODO WARN
|
|
return QPMS_SUCCESS;
|
|
}
|
|
|
|
/// Conversion from universal VSWF index u to type and the traditional layout index (\a l > 0).
|
|
/** Does *not* allow for longitudinal waves. */
|
|
static inline qpms_errno_t qpms_uvswfi2ty(qpms_uvswfi_t u,
|
|
qpms_vswf_type_t *t, qpms_y_t *y) {
|
|
*t = (qpms_vswf_type_t)(u & 3);
|
|
*y = u / 4 - 1;
|
|
if (QPMS_UNLIKELY(*t == 0 || *t == 3)) return QPMS_ERROR;
|
|
if (QPMS_UNLIKELY(*y < 0)) return QPMS_ERROR;
|
|
return QPMS_SUCCESS;
|
|
}
|
|
|
|
/// Conversion from universal VSWF index u to type and the traditional (vector) layout index.
|
|
/** *Does* allow for longitudinal waves, but returns an error if l == 0.
|
|
* (The only legit VSWF with l = 0 is the longitudinal one; u == 0.)
|
|
*/
|
|
static inline qpms_errno_t qpms_uvswfi2ty_l(qpms_uvswfi_t u,
|
|
qpms_vswf_type_t *t, qpms_y_t *y) {
|
|
*t = (qpms_vswf_type_t)(u & 3);
|
|
*y = u / 4 - 1;
|
|
if (QPMS_UNLIKELY(*t == 3)) return QPMS_ERROR;
|
|
if (QPMS_UNLIKELY(*y < 0)) return QPMS_ERROR;
|
|
return QPMS_SUCCESS;
|
|
}
|
|
|
|
/// Extract degree \a m from an universal VSWF index \a u.
|
|
static inline qpms_m_t qpms_uvswfi2m(qpms_uvswfi_t u) {
|
|
qpms_vswf_type_t t; qpms_m_t m; qpms_l_t n;
|
|
qpms_uvswfi2tmn(u, &t,&m,&n);
|
|
return m;
|
|
}
|
|
|
|
|
|
#ifdef __cplusplus
|
|
}
|
|
#endif
|
|
#endif //QPMS_INDEXING_H
|