From 3d83d174bdd744587bb133f7efd8aa2918686ac2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 12 Sep 2018 20:12:49 +0300 Subject: [PATCH] Fix points2d_rordered_frompoints last layer r; dump lattices intests/ewalds.c Former-commit-id: 1addca6a1c240310e86213c283df4162b2536b87 --- qpms/lattices2d.c | 10 ++++++++-- tests/ewalds.c | 38 +++++++++++++++++++++++++++++++++++--- 2 files changed, 43 insertions(+), 5 deletions(-) diff --git a/qpms/lattices2d.c b/qpms/lattices2d.c index fd0f4bd..9a2bbe2 100644 --- a/qpms/lattices2d.c +++ b/qpms/lattices2d.c @@ -66,9 +66,14 @@ points2d_rordered_t points2d_rordered_annulus(const points2d_rordered_t *orig, ptrdiff_t imin, imax; imin = points2d_rordered_locate_r(orig, minr); imax = points2d_rordered_locate_r(orig, maxr); + // TODO check + if(imax >= orig->nrs) --imax; + if(imax < 0) goto nothing; + // END TODO if (!inc_minr && (orig->rs[imin] <= minr)) ++imin; if (!inc_maxr && (orig->rs[imax] >= maxr)) --imax; if (imax < imin) { // it's empty +nothing: p.nrs = 0; p.base = NULL; p.rs = NULL; @@ -148,7 +153,7 @@ static points2d_rordered_t *points2d_rordered_frompoints_c(point2d *orig_base, c for (size_t i = 0; i < nmemb; ++i) { rcur = prn(p->base[i]); if (rcur > rcurmax) { - p->rs[ri] = rcursum / rcurcount; // average of the accrued r's within tolerance + p->rs[ri] = rcursum / (double) rcurcount; // average of the accrued r's within tolerance ++ri; p->r_offsets[ri] = i; //r_offsets[ri-1] + rcurcount (is the same) rcurcount = 0; @@ -158,6 +163,7 @@ static points2d_rordered_t *points2d_rordered_frompoints_c(point2d *orig_base, c rcursum += rcur; ++rcurcount; } + p->rs[ri] = rcursum / (double) rcurcount; p->r_offsets[rcount] = nmemb; return p; @@ -178,7 +184,7 @@ points2d_rordered_t *points2d_rordered_shift(const points2d_rordered_t *orig, co point2d * shifted = malloc(n * sizeof(point2d)); for(size_t i = 0; i < n; ++i) - shifted[i] = cart2_add(orig->base[i], shift); + shifted[i] = cart2_add(orig->base[i+orig->r_offsets[0]], shift); return points2d_rordered_frompoints_c(shifted, n,rtol, atol, false); } diff --git a/tests/ewalds.c b/tests/ewalds.c index 09e91c7..a3e2c00 100644 --- a/tests/ewalds.c +++ b/tests/ewalds.c @@ -37,7 +37,8 @@ typedef struct ewaldtest_triang_results { ewaldtest_triang_params paramslist[] = { - { 3, {1.1, 0.23}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_HORIZONTAL}, +// lMax, beta, k, a, eta, maxR, maxK, csphase, orientation +#if 0 { 3, {1.1, 0.23}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL}, { 3, {1.1, 0.23}, 2.3, 0.97, 0.5, 30, 30, -1., TRIANGULAR_VERTICAL}, { 3, {1.1, 0.23}, 2.3, 0.97, 0.9, 30, 30, 1., TRIANGULAR_VERTICAL}, @@ -50,6 +51,17 @@ ewaldtest_triang_params paramslist[] = { { 6, {1.1, 0.23}, 2.3, 0.97, 4.5, 40, 40, 1., TRIANGULAR_VERTICAL}, { 6, {1.1, 0.23}, 2.3, 0.97, 2.3, 100, 100, 1., TRIANGULAR_VERTICAL}, { 6, {1.1, 0.23}, 2.3, 0.97, 2.9, 100, 100, 1., TRIANGULAR_VERTICAL}, +#endif + { 2, {1.1, 0.23}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL}, + { 2, {-1.1, -0.23}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL}, + { 2, {0, 1.1}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL}, + { 2, {0, -1.1}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL}, + { 2, {-1.1, 0}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 0}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 0}, 2.3, 0.97, 0.5, 10, 10, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 0}, 2.3, 0.97, 0.5, 5, 5, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 0}, 2.3, 0.97, 0.5, 2, 2, 1., TRIANGULAR_VERTICAL}, + // end: // { 0, {0, 0}, 0, 0, 0, 0, 0, 0, 0} }; @@ -65,6 +77,18 @@ void ewaldtest_triang_results_free(ewaldtest_triang_results *r) { free(r); } + +void dump_points2d_rordered(const points2d_rordered_t *ps, char *filename) { + FILE *f = fopen(filename, "w"); + for (size_t i = 0; i < ps->nrs; ++i) { + fprintf(f, "# r = %.16g\n", ps->rs[i]); + for (ptrdiff_t j = ps->r_offsets[i]; j < ps->r_offsets[i+1]; ++j) + fprintf(f, "%.16g %.16g\n", ps->base[j].x, ps->base[j].y); + } + fclose(f); +} + + ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p); int main() { @@ -74,8 +98,8 @@ int main() { ewaldtest_triang_results *r = ewaldtest_triang(p); // TODO print per-test header here printf("===============================\n"); - printf("Kmax = %g, Rmax = %g, lMax = %d, eta = %g, k = %g, beta = (%g,%g), csphase = %g\n", - p.maxK, p.maxR, p.lMax, p.eta, p.k, p.beta.x, p.beta.y, p.csphase); + printf("a = %g, K = %g, Kmax = %g, Rmax = %g, lMax = %d, eta = %g, k = %g, beta = (%g,%g), csphase = %g\n", + p.a, 4*M_PI/sqrt(3)/p.a, p.maxK, p.maxR, p.lMax, p.eta, p.k, p.beta.x, p.beta.y, p.csphase); printf("sigma0: %.16g%+.16gj\n", creal(r->sigma0), cimag(r->sigma0)); for (qpms_l_t n = 0; n <= p.lMax; ++n) { for (qpms_m_t m = -n; m <= n; ++m){ @@ -101,6 +125,9 @@ int main() { } +int ewaldtest_counter = 0; + + ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) { const double a = p.a; //const double a = p.h * sqrt(3); @@ -134,6 +161,10 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) { points2d_rordered_t *Kpoints_plus_beta = points2d_rordered_shift(&(Klg->ps), p.beta, 8*DBL_EPSILON, 8*DBL_EPSILON); + char filename[BUFSIZ]; + sprintf(filename, "betalattice_%d.out", ewaldtest_counter); + dump_points2d_rordered(Kpoints_plus_beta, filename); + point2d particle_shift = {0,0}; // TODO make this a parameter if (0!=ewald32_sigma_long_shiftedpoints(results->sigmas_long, @@ -192,5 +223,6 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) { qpms_ewald32_constants_free(c); triangular_lattice_gen_free(Klg); triangular_lattice_gen_free(Rlg); + ++ewaldtest_counter; return results; }