ewald.h

Go to the documentation of this file.

 
#ifndef EWALD_H
#define EWALD_H
#include <gsl/gsl_sf_result.h>
#include <stdlib.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_errno.h>
#include <math.h> // for inlined lilgamma
#include <complex.h>
#include "qpms_types.h"
#include "lattices.h"
 
 
typedef enum {
    QPMS_EWALD_LONG_RANGE = 1,
    QPMS_EWALD_SHORT_RANGE = 2,
    QPMS_EWALD_0TERM = 4,
    QPMS_EWALD_FULL = QPMS_EWALD_LONG_RANGE | QPMS_EWALD_SHORT_RANGE | QPMS_EWALD_0TERM,
} qpms_ewald_part;
 
 
gsl_error_handler_t IgnoreUnderflowsGSLErrorHandler;
 
 
 
typedef struct qpms_ewald3_constants_t { 
    qpms_l_t lMax;
    qpms_y_t nelem_sc;
    qpms_l_t *s1_jMaxes;
    complex double **s1_constfacs; // indices [y][j] where j is same as in [1, (4.5)]
    /* These are the actual numbers now: (in the EWALD32_CONSTANTS_AGNOSTIC version)
     * for m + n EVEN:
     *
     * s1_constfacs[y(m,n)][j] = 
     *
     *   -2 * I**(n+1) * sqrt(π) * ((n-m)/2)! * ((n+m)/2)! * (-1)**j 
     *   -----------------------------------------------------------
     *              j! * ((n-m)/2 - j)! * ((n+m)/2 + j)!
     *
     * for m + n ODD:
     *
     * s1_constfacs[y(m,n)][j] = 0
     */
    complex double *s1_constfacs_base; 
    // similarly for the 1D z-axis aligned case; now the indices are [n][j] (as m == 0)
 
    
    // TODO indexing mechanisms
    
    complex double **S1_constfacs; // indices [y][j] where j is same as in [1, (4.5)]
    /* These are the actual numbers now: (in the EWALD32_CONSTANTS_AGNOSTIC version)
     * for m + n EVEN:
     *
     * S1_constfacs[y(m,n)][x(j,s)] = 
     *
     *   -2 * I**(n+1) * sqrt(π) * ((n-m)/2)! * ((n+m)/2)! * (-1)**j /   j    \
     *   ----------------------------------------------------------- |        |
     *              j! * ((n - m - s)/2)! * ((n + m - s)/2)!         \ 2j - s /
     *
     * for m + n ODD:
     *
     * S1_constfacs[y(m,n)][j] = 0
     */
    complex double *S1_constfacs_base; 
 
    complex double **s1_constfacs_1Dz; 
    /* These are the actual numbers now:
     * s1_constfacs_1Dz[n][j] =
     *   
     *     -I**(n+1) (-1)**j * n!
     *   --------------------------
     *   j! * 2**(2*j) * (n - 2*j)!
     */
    complex double *s1_constfacs_1Dz_base; 
 
    double *legendre0; /* now with GSL_SF_LEGENDRE_NONE normalisation, because this is what is
                * what the multipliers from translations.c count with.
                */
    double *legendre_plus1; // needed? TODO; in any case, nonzero only for m=0
    double *legendre_minus1; // needed? TODO; in any case, nonzero only for m=0
    gsl_sf_legendre_t legendre_normconv;
    int legendre_csphase;       /* 1 or -1; csphase of the Legendre polynomials saved in legendre0 etc.
                    This is because I dont't actually consider this fixed in
                    translations.c */
 
} qpms_ewald3_constants_t;
 
qpms_ewald3_constants_t *qpms_ewald3_constants_init(qpms_l_t lMax, int csphase);
void qpms_ewald3_constants_free(qpms_ewald3_constants_t *);
 
 
 
typedef struct qpms_csf_result { 
  complex double val; 
  double err; 
} qpms_csf_result;
 
 
// [1, (A.9)]
static inline complex double lilgamma(double t) {
  t = fabs(t);
  if (t >= 1) 
    return sqrt(t*t - 1);
  else
    return -I * sqrt(1 - t*t);
}
 
// [1, (A.8)], complex version of lilgamma()
static inline complex double clilgamma(complex double z) {
    complex double a1 = z - 1, a2 = z + 1;
    // ensure  -pi/2 < arg(z + 1) < 3*pi/2
    if (creal(a2) < 0 && cimag(a2) <= 0) 
        a2 = -csqrt(a2);
    else 
        a2 = csqrt(a2);
    // ensure -3*pi/2 < arg(z - 1) < pi/2
    if (creal(a1) < 0 && cimag(a1) >= 0)
        a1 = -csqrt(a1);
    else
        a1 = csqrt(a1);
    return a1 * a2;
}
 
 
int cx_gamma_inc_series_e(double a, complex double z, qpms_csf_result * result);
 
 
int cx_gamma_inc_CF_e(double a, complex double z, qpms_csf_result * result);
 
 
int complex_gamma_inc_e(double a, complex double x,
 
    int m,
    qpms_csf_result *result);
 
 
int complex_expint_n_e(int n, complex double x, qpms_csf_result *result);
 
 
int hyperg_2F2_series(double a, double b, double c, double d, double x, 
        gsl_sf_result *result);
 
#if 0
// The integral from (4.6); maybe should be static and not here.
int ewald32_sr_integral(double r, double k, double n, double eta, double *result, double *err, gsl_integration_workspace *workspace);
#endif
 
 
void ewald3_2_sigma_long_Delta(complex double *target, double *target_err, int maxn, complex double x,
           int xbranch, complex double z);
 
 
void ewald3_2_sigma_long_Delta_recurrent(complex double *target, double *target_err, int maxn, complex double x,
           int xbranch, complex double z, _Bool bigimz);
 
 
void ewald3_2_sigma_long_Delta_series(complex double *target, double *target_err, int maxn, complex double x,
           int xbranch, complex double z);
 
// General functions acc. to [2], sec. 4.6 – currently valid for 2D and 1D lattices in 3D space
 
int ewald3_sigma0(complex double *result, 
            double *err, 
        const qpms_ewald3_constants_t *c, 
        double eta, 
        complex double wavenumber 
);
 
int ewald3_sigma_short(
        complex double *target_sigmasr_y, 
        double *target_sigmasr_y_err, 
        const qpms_ewald3_constants_t *c, 
        double eta, 
            complex double wavenumber, 
 
        LatticeDimensionality latdim,
 
        PGen *pgen_R,
 
        bool pgen_generates_shifted_points,
        cart3_t k,
        cart3_t particle_shift
        );
 
int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, depending on latdim
        complex double *target_sigmalr_y, 
        double *target_sigmalr_y_err, 
        const qpms_ewald3_constants_t *c, 
        double eta, 
            complex double wavenumber, 
            double unitcell_volume, 
        LatticeDimensionality latdim,
 
        PGen *pgen_K, 
 
        bool pgen_generates_shifted_points,
        cart3_t k,
        cart3_t particle_shift
        );
        
#endif //EWALD_H