qpms/normalisation.h

   1 /*! \file normalisation.h
   2  * \brief Convention-dependent coefficients for VSWFs.
   3  *
   4  * See also @ref qpms_normalisation_t and @ref vswf_conventions.
   5  */
   6 #ifndef NORMALISATION_H
   7 #define NORMALISATION_H
   8
   9 #include "qpms_types.h"
  10 #include "qpms_error.h"
  11 #include <math.h>
  12 #include <complex.h>
  13 #include "indexing.h"
  14 #include "optim.h"
  15
  16 /// Returns the (real positive) common norm factor of a given normalisation compared to the reference convention.
  17 /** Does NOT perform the inversion if QPMS_NORMALISATION_INVERSE is set. */
  18 static inline double qpms_normalisation_normfactor(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
  19         switch (QPMS_EXPECT(norm & QPMS_NORMALISATION_NORM_BITS, norm & QPMS_NORMALISATION_DEFAULT)) {
  20                 case QPMS_NORMALISATION_NORM_POWER:
  21                         return 1;
  22                 case QPMS_NORMALISATION_NORM_SPHARM:
  23                         return sqrt(l*(l+1));
  24                 case QPMS_NORMALISATION_NORM_NONE: // TODO more precision
  25                         return sqrt(l*(l+1) * 4*M_PI / (2*l+1)) *
  26                                 exp(0.5*(lgamma(l+m+1) - lgamma(l-m+1)));
  27                 default:
  28                         QPMS_WTF;
  29         }
  30 }
  31
  32
  33
  34 /// Returns the factors of a magnetic basis VSWF of a given convention compared to the reference convention.
  35 /**
  36  * This version ignores the Condon-Shortley phase bit (perhaps because the Condon-Shortley
  37  * phase is already taken into account in a `gsl_sf_legendre_*_e()` call.)
  38  */
  39 static inline complex double qpms_normalisation_factor_M_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
  40         complex double fac = qpms_normalisation_normfactor(norm, l, m);
  41         if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_M_MINUS)) fac *= -1;
  42         if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_M_I)) fac *= I;
  43         if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac;
  44         return fac;
  45 }
  46
  47
  48 /// Returns the factors of a magnetic basis VSWF of a given convention compared to the reference convention.
  49 /**
  50  * This version takes into account the Condon-Shortley phase bit.
  51  * Do not use if the C.-S. has already been taken into account e.g. in
  52  * a `gsl_sf_legendre_*_e()` call.
  53  */
  54 static inline complex double qpms_normalisation_factor_M(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
  55         complex double fac = qpms_normalisation_factor_M_noCS(norm, l, m);
  56         return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac;
  57 }
  58
  59
  60 /// Returns the factors of a electric basis VSWF of a given convention compared to the reference convention.
  61 /**
  62  * This version ignores the Condon-Shortley phase bit (perhaps because the Condon-Shortley
  63  * phase is already taken into account in a `gsl_sf_legendre_*_e()` call.)
  64  */
  65 static inline complex double qpms_normalisation_factor_N_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
  66         complex double fac = qpms_normalisation_normfactor(norm, l, m);
  67         if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_N_MINUS)) fac *= -1;
  68         if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_N_I)) fac *= I;
  69         if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac;
  70         return fac;
  71 }
  72
  73
  74 /// Returns the factors of a electric basis VSWF of a given convention compared to the reference convention.
  75 /**
  76  * This version takes into account the Condon-Shortley phase bit.
  77  * Do not use if the C.-S. has already been taken into account e.g. in
  78  * a `gsl_sf_legendre_*_e()` call.
  79  */
  80 static inline complex double qpms_normalisation_factor_N(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
  81         complex double fac = qpms_normalisation_factor_N_noCS(norm, l, m);
  82         return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac;
  83 }
  84
  85
  86 /// Returns the factors of a electric basis VSWF divided by the factor of a magnetic VWFS of a given convention, compared to the reference one.
  87 static inline complex double qpms_normalisation_factor_N_M(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
  88         return qpms_normalisation_factor_N_noCS(norm, l, m)
  89                 / qpms_normalisation_factor_M_noCS(norm, l, m);
  90 }
  91
  92
  93 /// Returns the factors of a longitudinal basis VSWF of a given convention compared to the reference convention.
  94 /**
  95  * This version ignores the Condon-Shortley phase bit (perhaps because the Condon-Shortley
  96  * phase is already taken into account in a `gsl_sf_legendre_*_e()` call.)
  97  */
  98 static inline complex double qpms_normalisation_factor_L_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
  99         complex double fac = qpms_normalisation_normfactor(norm, l, m);
 100         if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_L_MINUS)) fac *= -1;
 101         if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_L_I)) fac *= I;
 102         if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac;
 103         return fac;
 104 }
 105
 106 /// Returns the factors of a longitudinal basis VSWF of a given convention compared to the reference convention.
 107 /**
 108  * This version takes into account the Condon-Shortley phase bit.
 109  * Do not use if the C.-S. has already been taken into account e.g. in
 110  * a `gsl_sf_legendre_*_e()` call.
 111  */
 112 static inline complex double qpms_normalisation_factor_L(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
 113         complex double fac = qpms_normalisation_factor_L_noCS(norm, l, m);
 114         return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac;
 115 }
 116
 117 /// Returns the factors of a basis VSWF of a given convention compared to the reference convention.
 118 static inline complex double qpms_normalisation_factor_uvswfi(const qpms_normalisation_t norm, qpms_uvswfi_t ui) {
 119         qpms_vswf_type_t t; qpms_m_t m; qpms_l_t l;
 120         qpms_uvswfi2tmn(ui, &t, &m, &l);
 121         switch(t) {
 122                 case QPMS_VSWF_MAGNETIC:
 123                         return qpms_normalisation_factor_M(norm, l, m);
 124                 case QPMS_VSWF_ELECTRIC:
 125                         return qpms_normalisation_factor_N(norm, l, m);
 126                 case QPMS_VSWF_LONGITUDINAL:
 127                         return qpms_normalisation_factor_L(norm, l, m);
 128                 default:
 129                         QPMS_WTF;
 130         }
 131 }
 132
 133
 134 /// Returns normalisation flags corresponding to the dual spherical harmonics / waves.
 135 /**
 136  * This reverses the normalisation factors returned by qpms_normalisation_factor_*
 137  * and conjugates the asimuthal part for complex spherical harmonics,
 138  * \f$ e^{\pm im\phi} \leftrightarrow e^{\mp im\phi} \f$.
 139  */
 140 static inline qpms_normalisation_t qpms_normalisation_dual(qpms_normalisation_t norm) {
 141         norm ^= QPMS_NORMALISATION_INVERSE;
 142         if (!(norm & QPMS_NORMALISATION_SPHARM_REAL))
 143                 norm ^= QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE;
 144         return norm;
 145 }
 146
 147 /// Returns the asimuthal part of a spherical harmonic.
 148 /** Returns \f[ e^{im\phi} \f] for standard complex spherical harmonics,
 149  * \f[ e^{-im\phi \f] for complex spherical harmonics
 150  * and QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE set.
 151  *
 152  * For real spherical harmonics, this gives
 153  * \f[
 154  *      \sqrt{2}\cos{m \phi} \quad \mbox{if } m>0, \\
 155  *      \sqrt{2}\sin{m \phi} \quad \mbox{if } m<0, \\
 156  *      0 \quad \mbox{if } m>0. \\
 157  * \f]
 158  */
 159 static inline complex double qpms_spharm_azimuthal_part(qpms_normalisation_t norm, qpms_m_t m, double phi) {
 160         switch(QPMS_EXPECT(norm, QPMS_NORMALISATION_DEFAULT)
 161                         & (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) {
 162                 case 0:
 163                         return cexp(I*m*phi);
 164                 case QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE:
 165                         return cexp(-I*m*phi);
 166                 case QPMS_NORMALISATION_SPHARM_REAL:
 167                         if (m > 0) return M_SQRT2 * cos(m*phi);
 168                         else if (m < 0) return M_SQRT2 * sin(m*phi);
 169                         else return 1.;
 170                 default:
 171                         QPMS_WTF;
 172         }
 173 }
 174
 175 /// Returns derivative of the asimuthal part of a spherical harmonic divided by \a m.
 176 /**
 177  *
 178  * This is used to evaluate the VSWFs together with the \a pi member array of the
 179  * qpms_pitau_t structure.
 180  *
 181  * Returns \f[ i e^{im\phi} \f] for standard complex spherical harmonics,
 182  * \f[-i e^{-i\phi \f] for complex spherical harmonics
 183  * and QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE set.
 184  *
 185  * For real spherical harmonics, this gives
 186  * \f[
 187  *      -\sqrt{2}\sin{m \phi} \quad \mbox{if } m>0, \\
 188  *      \sqrt{2}\cos{m \phi} \quad \mbox{if } m<0, \\
 189  *      -1 \quad \mbox{if } \mbox{if }m=0. \\
 190  * \f]
 191  *
 192  * (The value returned for \f$ m = 0 \f$ should not actually be used for
 193  * anything except for multiplying by zero.)
 194  *
 195  *
 196  */
 197 static inline complex double qpms_spharm_azimuthal_part_derivative_div_m(qpms_normalisation_t norm, qpms_m_t m, double phi) {
 198         if(m==0) return 0;
 199         switch(QPMS_EXPECT(norm, QPMS_NORMALISATION_DEFAULT)
 200                         & (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) {
 201                 case 0:
 202                         return I*cexp(I*m*phi);
 203                 case QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE:
 204                         return -I*cexp(-I*m*phi);
 205                 case QPMS_NORMALISATION_SPHARM_REAL:
 206                         if (m > 0) return -M_SQRT2 * sin(m*phi);
 207                         else if (m < 0) return M_SQRT2 * cos(m*phi);
 208                         else return -1;
 209                 default:
 210                         QPMS_WTF;
 211         }
 212 }
 213
 214 #if 0 // legacy code moved from qpms_types.h. TODO cleanup
 215 /// Returns the normalisation convention code without the Condon-Shortley phase.
 216 static inline int qpms_normalisation_t_normonly(qpms_normalisation_t norm) {
 217         return norm & (~QPMS_NORMALISATION_T_CSBIT);
 218 }
 219
 220 /* Normalisation of the spherical waves is now scattered in at least three different files:
 221  * here, we have the norm in terms of radiated power of outgoing wave.
 222  * In file legendre.c, function qpms_pitau_get determines the norm used in the vswf.c
 223  * spherical vector wave norms. The "dual" waves in vswf.c use the ..._abssquare function below.
 224  * In file translations.c, the normalisations are again set by hand using the normfac and lognormfac
 225  * functions.
 226  */
 227 #include <math.h>
 228 #include <assert.h>
 229 // relative to QPMS_NORMALISATION_KRISTENSSON_CS, i.e.
 230 // P_l^m[normtype] = P_l^m[Kristensson]
 231 static inline double qpms_normalisation_t_factor(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
 232         int csphase = qpms_normalisation_t_csphase(norm);
 233         norm = qpms_normalisation_t_normonly(norm);
 234         double factor;
 235         switch (norm) {
 236                 case QPMS_NORMALISATION_KRISTENSSON:
 237                         factor = 1.;
 238                         break;
 239                 case QPMS_NORMALISATION_TAYLOR:
 240                         factor = sqrt(l*(l+1));
 241                         break;
 242                 case QPMS_NORMALISATION_NONE:
 243                         factor = sqrt(l*(l+1) * 4 * M_PI / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1)));
 244                         break;
 245 #ifdef USE_XU_ANTINORMALISATION // broken probably in legendre.c
 246                 case QPMS_NORMALISATION_XU:
 247                         factor = sqrt(4 * M_PI) / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1));
 248                         break;
 249 #endif
 250                 default:
 251                         assert(0);
 252         }
 253         factor *= (m%2)?(-csphase):1;
 254         return factor;
 255 }
 256
 257
 258 // TODO move elsewhere
 259 static inline double qpms_normalisation_t_factor_abssquare(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
 260         norm = qpms_normalisation_t_normonly(norm);
 261         switch (norm) {
 262                 case QPMS_NORMALISATION_KRISTENSSON:
 263                         return 1.;
 264                         break;
 265                 case QPMS_NORMALISATION_TAYLOR:
 266                         return l*(l+1);
 267                         break;
 268                 case QPMS_NORMALISATION_NONE:
 269                         return l*(l+1) * 4 * M_PI / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1));
 270                         break;
 271 #ifdef USE_XU_ANTINORMALISATION // broken probably in legendre.c
 272                 case QPMS_NORMALISATION_XU:
 273                         {
 274                           double fac = sqrt(4 * M_PI) / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1));
 275                           return fac * fac;
 276                         }
 277                         break;
 278 #endif
 279                 default:
 280                         assert(0);
 281                         return NAN;
 282         }
 283 }
 284 #endif
 285
 286 #endif //NORMALISATION_H