qpms/translations.h

   1 /*! \file translations.h
   2  *  \brief VSWF translation operator.
   3  *
   4  * ### Argument conventions
   5  *
   6  * A single wave with indices mu, nu is re-expanded at kdlj into waves with indices m, n,
   7  * i.e. in the following functions, the first arguments over which one sums (multiplied
   8  * by the waves with new origin).
   9  *
  10  * HOWEVER, this means that if a field has an expansion with coeffs a(mu, nu)
  11  * at the original origin, with the basis at the new origin, the coeffs will be
  12  * a(m, n) = \sum_{mu,nu} A(m, n, mu, nu) a(mu, nu).
  13  *
  14  * With qpms_trans_calculator_get_AB_arrays_buf (and other functions from *AB_arrays*
  15  * family), one can choose the stride. And it seems that the former stride argument (now called
  16  * destride) and the latter (now called srcstride) are connected to (m,n) and (mu,nu) indices,
  17  * respectively. Seems consistent.
  18  *
  19  *
  20  * #### r_ge_d argument:
  21  *
  22  * If r_ge_d == true, the translation coefficients are calculated using regular bessel functions,
  23  * regardless of what J argument is.
  24  *
  25  *
  26  */
  27 #ifndef QPMS_TRANSLATIONS_H
  28 #define QPMS_TRANSLATIONS_H
  29 #include "vectors.h"
  30 #include "qpms_types.h"
  31 #include <complex.h>
  32 #include <stdbool.h>
  33 #include <stddef.h>
  34
  35 #if defined LATTICESUMS32 || defined LATTICESUMS31
  36 #include "ewald.h"
  37 #endif
  38
  39
  40 complex double qpms_trans_single_A(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
  41                 bool r_ge_d, qpms_bessel_t J);
  42
  43 complex double qpms_trans_single_B(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
  44                 bool r_ge_d, qpms_bessel_t J);
  45
  46 /// Structure holding the constant factors in normalisation operators.
  47 /**
  48  * The VSWF translation operator elements are rather complicated linear
  49  * combinations of Bessel functions and spherical harmonics.
  50  * The preceding factors are rather complicated but need to be calculated
  51  * (for a single normalisation convention)
  52  * only once and then recycled during the operator evaluation for different
  53  * translations.
  54  *
  55  * This structure is initialised with qpms_trans_calculator_t() and
  56  * holds all the constant factors up to a truncation
  57  * degree \a lMax.
  58  *
  59  * The destructor function is qpms_trans_calculator_free().
  60  *
  61  * If Ewald sums are enabled at build, it also holds the constant
  62  * factors useful for lattice sums of translation operator.
  63  */
  64 typedef struct qpms_trans_calculator {
  65         qpms_normalisation_t normalisation;
  66         qpms_l_t lMax;
  67         qpms_y_t nelem;
  68         complex double **A_multipliers;
  69         complex double **B_multipliers;
  70 #if 0
  71         // Normalised values of the Legendre functions and derivatives
  72         // for θ == π/2, i.e. for the 2D case.
  73         double *leg0;
  74         double *pi0;
  75         double *tau0;
  76         // Spherical Bessel function coefficients:
  77         // TODO
  78 #endif
  79
  80 #if defined LATTICESUMS32 || defined LATTICESUMS31
  81         qpms_ewald3_constants_t *e3c;
  82 #endif
  83         double *legendre0; // Zero-argument Legendre functions – this might go outside #ifdef in the end...
  84 } qpms_trans_calculator;
  85
  86 /// Initialise a qpms_trans_calculator_t instance for a given convention and truncation degree.
  87 qpms_trans_calculator *qpms_trans_calculator_init(qpms_l_t lMax, ///< Truncation degree.
  88                 qpms_normalisation_t nt
  89                 );
  90 /// Destructor for qpms_trans_calculator_t.
  91 void qpms_trans_calculator_free(qpms_trans_calculator *);
  92
  93 complex double qpms_trans_calculator_get_A(const qpms_trans_calculator *c,
  94                 qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
  95                 bool r_ge_d, qpms_bessel_t J);
  96 complex double qpms_trans_calculator_get_B(const qpms_trans_calculator *c,
  97                 qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
  98                 bool r_ge_d, qpms_bessel_t J);
  99 int qpms_trans_calculator_get_AB_p(const qpms_trans_calculator *c,
 100                 complex double *Adest, complex double *Bdest,
 101                 qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
 102                 bool r_ge_d, qpms_bessel_t J);
 103 int qpms_trans_calculator_get_AB_arrays(const qpms_trans_calculator *c,
 104                 complex double *Adest, complex double *Bdest,
 105                 size_t deststride, size_t srcstride,
 106                 csph_t kdlj, bool r_ge_d, qpms_bessel_t J);
 107
 108 // TODO update the types later
 109 complex double qpms_trans_calculator_get_A_ext(const qpms_trans_calculator *c,
 110                 int m, int n, int mu, int nu, complex double kdlj_r,
 111                 double kdlj_th, double kdlj_phi, int r_ge_d, int J);
 112
 113 complex double qpms_trans_calculator_get_B_ext(const qpms_trans_calculator *c,
 114                 int m, int n, int mu, int nu, complex double kdlj_r,
 115                 double kdlj_th, double kdlj_phi, int r_ge_d, int J);
 116
 117 int qpms_trans_calculator_get_AB_p_ext(const qpms_trans_calculator *c,
 118                 complex double *Adest, complex double *Bdest,
 119                 int m, int n, int mu, int nu, complex double kdlj_r,
 120                 double kdlj_th, double kdlj_phi, int r_ge_d, int J);
 121
 122 int qpms_trans_calculator_get_AB_arrays_ext(const qpms_trans_calculator *c,
 123                 complex double *Adest, complex double *Bdest,
 124                 size_t deststride, size_t srcstride,
 125                 complex double kdlj_r, double kdlj_theta, double kdlj_phi,
 126                 int r_ge_d, int J);
 127
 128 // Convenience functions using VSWF base specs
 129 qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator *c,
 130                 complex double *target,
 131                 /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
 132                 const qpms_vswf_set_spec_t *destspec, size_t deststride,
 133                 /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
 134                 const qpms_vswf_set_spec_t *srcspec, size_t srcstride,
 135                 csph_t kdlj, bool r_ge_d, qpms_bessel_t J);
 136
 137 /// Version with \a k and cartesian particle positions
 138 /// and with automatic \a r_ge_d = `false`.
 139 qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p(
 140                 const qpms_trans_calculator *c,
 141                 complex double *target,
 142                 /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
 143                 const qpms_vswf_set_spec_t *destspec, size_t deststride,
 144                 /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
 145                 const qpms_vswf_set_spec_t *srcspec, size_t srcstride,
 146                 complex double k, cart3_t destpos, cart3_t srcpos,
 147                 qpms_bessel_t J
 148                 /// Workspace has to be at least 2 * c->neleme**2 long
 149                 );
 150
 151 #ifdef LATTICESUMS32
 152 // for the time being there are only those relatively low-level quick&dirty functions
 153 // according to what has been implemented from ewald.h;
 154 // TODO more high-level functions with more advanced lattice generators etc. (after
 155 // the prerequisities from lattices2d.h are implememted)
 156
 157
 158 int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
 159                 complex double *Adest, double *Aerr,
 160                 complex double *Bdest, double *Berr,
 161                 const ptrdiff_t deststride, const ptrdiff_t srcstride,
 162                 const double eta, const complex double wavenumber,
 163                 cart2_t b1, cart2_t b2,
 164                 const cart2_t beta,
 165                 const cart3_t particle_shift,
 166                 double maxR, double maxK
 167                 );
 168
 169 int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
 170                 complex double *Adest, double *Aerr,
 171                 complex double *Bdest, double *Berr,
 172                 const ptrdiff_t deststride, const ptrdiff_t srcstride,
 173                 const double eta, const complex double wavenumber,
 174                 cart2_t b1, cart2_t b2,
 175                 const cart2_t beta,
 176                 const cart3_t particle_shift,
 177                 double maxR, double maxK,
 178                 qpms_ewald_part parts
 179                 );
 180
 181 // Convenience functions using VSWF base specs
 182 qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculator *c,
 183                 complex double *target, double *err,
 184                 /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
 185                 const qpms_vswf_set_spec_t *destspec, size_t deststride,
 186                 /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
 187                 const qpms_vswf_set_spec_t *srcspec, size_t srcstride,
 188                 const double eta, const complex double wavenumber,
 189                 cart2_t b1, cart2_t b2,
 190                 const cart2_t beta,
 191                 const cart3_t particle_shift,
 192                 double maxR, double maxK
 193                 );
 194
 195 qpms_errno_t qpms_trans_calculator_get_trans_array_e32_e(const qpms_trans_calculator *c,
 196                 complex double *target, double *err,
 197                 /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
 198                 const qpms_vswf_set_spec_t *destspec, size_t deststride,
 199                 /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
 200                 const qpms_vswf_set_spec_t *srcspec, size_t srcstride,
 201                 const double eta, const complex double wavenumber,
 202                 cart2_t b1, cart2_t b2,
 203                 const cart2_t beta,
 204                 const cart3_t particle_shift,
 205                 double maxR, double maxK,
 206                 qpms_ewald_part parts
 207                 );
 208
 209 #endif //LATTICESUMS32
 210
 211 #ifdef LATTICESUMS31
 212 // e31z means that the particles are positioned along the z-axis;
 213 // their positions and K-values are then denoted by a single z-coordinate
 214 int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_trans_calculator *c,
 215                 complex double *Adest, double *Aerr,
 216                 complex double *Bdest, double *Berr,
 217                 const ptrdiff_t deststride, const ptrdiff_t srcstride,
 218                 const double eta, const double k,
 219                 const double unitcell_area, // just lattice period
 220                 const size_t nRpoints, const cart2_t *Rpoints,
 221                 const size_t nKpoints, const cart2_t *Kpoints,
 222                 const double beta,
 223                 const double particle_shift
 224                 );
 225 #endif
 226
 227
 228
 229
 230 #ifdef QPMS_COMPILE_PYTHON_EXTENSIONS
 231 #include <Python.h>
 232 #include <numpy/npy_common.h>
 233 int qpms_cython_trans_calculator_get_AB_arrays_loop(
 234                 const qpms_trans_calculator *c, qpms_bessel_t J, const int resnd,
 235                 int daxis, int saxis,
 236                 char *A_data, const npy_intp *A_shape, const npy_intp *A_strides,
 237                 char *B_data, const npy_intp *B_shape, const npy_intp *B_strides,
 238                 const char *r_data, const npy_intp *r_shape, const npy_intp *r_strides,
 239                 const char *theta_data, const npy_intp *theta_shape, const npy_intp *theta_strides,
 240                 const char *phi_data, const npy_intp *phi_shape, const npy_intp *phi_strides,
 241                 const char *r_ge_d_data, const npy_intp *r_ge_d_shape, const npy_intp *r_ge_d_strides);
 242
 243
 244 #endif //QPMS_COMPILE_PYTHON_EXTENSIONS
 245
 246
 247 #endif // QPMS_TRANSLATIONS_H