From 12c92a6423b3dc3c1a148296611413c512335d0a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 19 Sep 2018 10:23:18 +0300 Subject: [PATCH] Fix singular term Former-commit-id: 33a368c7126d4a8dfa5bdebaae3e70badfd6c2f5 --- qpms/apps/hexlattice_ewald.c | 10 +++++++++- qpms/translations.c | 2 +- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/qpms/apps/hexlattice_ewald.c b/qpms/apps/hexlattice_ewald.c index b944c0e..9dbe51d 100644 --- a/qpms/apps/hexlattice_ewald.c +++ b/qpms/apps/hexlattice_ewald.c @@ -73,13 +73,21 @@ int main (int argc, char **argv) { const double unitcell_area = s3*a*a/2.; const double rec_a = 4*M_PI/s3/a; + //fprintf(stderr, "%.16g\n", rec_a); + triangular_lattice_gen_t *Rlg = triangular_lattice_gen_init(a, rs_orientation, true, 0); triangular_lattice_gen_extend_to_steps(Rlg, RLAYERS); triangular_lattice_gen_t *Klg = triangular_lattice_gen_init(rec_a, reverseTriangularLatticeOrientation(rs_orientation), true, 0); triangular_lattice_gen_extend_to_steps(Klg, KLAYERS); + points2d_rordered_t Rselwo0 = points2d_rordered_annulus(&(Rlg->ps), 0, false, INFINITY, false); + const point2d *Rpoints = Rlg->ps.base; size_t nR = Rlg->ps.r_offsets[Rlg->ps.nrs]; + // TODO automatic skip of the zero point directly in translations.c + const point2d *Rpoints_wo0 = Rselwo0.base + Rselwo0.r_offsets[0]; + size_t nR_wo0 = Rselwo0.r_offsets[Rselwo0.nrs] - Rselwo0.r_offsets[0]; + assert(nR - nR_wo0 == 1); const point2d *Kpoints = Klg->ps.base; size_t nK = Klg->ps.r_offsets[Klg->ps.nrs]; @@ -131,7 +139,7 @@ int main (int argc, char **argv) { &(W[0][0][0][0][0]), err ? &(Werr[0][0][0][0][0]) : NULL, // Adest, Aerr, &(W[0][0][1][0][0]), err ? &(Werr[0][0][1][0][0]) : NULL, // Bdest, Berr, deststride, srcstride, - eta, k0_eff, unitcell_area, nR, Rpoints, nK, Kpoints, beta, pshift0 + eta, k0_eff, unitcell_area, nR_wo0, Rpoints_wo0, nK, Kpoints, beta, pshift0 ); // B<-B for(qpms_y_t desty = 0; desty < c->nelem; ++desty) diff --git a/qpms/translations.c b/qpms/translations.c index 564f4ad..69f52b7 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -1163,7 +1163,7 @@ int qpms_trans_calculator_get_AB_arrays_e32_both_points_and_shift(const qpms_tra const ptrdiff_t deststride, const ptrdiff_t srcstride, /* qpms_bessel_t J*/ // assume QPMS_HANKEL_PLUS const double eta, const double k, const double unitcell_area, - const size_t nRpoints, const cart2_t *Rpoints, + const size_t nRpoints, const cart2_t *Rpoints, // n.b. can't contain 0; TODO automatic recognition and skip const size_t nKpoints, const cart2_t *Kpoints, const cart2_t beta, const cart2_t particle_shift