From 60e1490d74b968959f5f632ffa03d1daa2478283 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sun, 26 Aug 2018 16:30:07 +0000 Subject: [PATCH] "Annulus view" for points2d_rordered_t; +prototypes for functions inspired by lattices2d.py (not implemented) Former-commit-id: da46b26573e3fd23b3cb3e8a0fe20c6bc3dff3e6 --- qpms/lattices.h | 108 ++++++++++++++++++++++++++ qpms/lattices2d.c | 36 +++++++++ tests/lattice/2d_triangular_lattice.c | 12 ++- 3 files changed, 155 insertions(+), 1 deletion(-) diff --git a/qpms/lattices.h b/qpms/lattices.h index df0d053..987655a 100644 --- a/qpms/lattices.h +++ b/qpms/lattices.h @@ -1,5 +1,17 @@ #ifndef LATTICES_H #define LATTICES_H + + +/* IMPORTANT TODO + * ============== + * + * The current content of this part (and the implementation) is extremely ugly. + * When I have some time, I have to rewrite this in the style of lattices2d.py + * + */ + + + #include #include #include @@ -20,6 +32,87 @@ static inline point2d point2d_fromxy(const double x, const double y) { return p; } +/* + * THE NICE PART (adaptation of lattices2d.py) + * =========================================== + * + * all the functions are prefixed with l2d_ + * convention for argument order: inputs, *outputs, add. params + */ + +#define BASIS_RTOL 1e-13 + +typedef enum { + OBLIQUE = 1, + RECTANGULAR = 2, + SQUARE = 4, + RHOMBIC = 5, + EQUILATERAL_TRIANGULAR = 3, + RIGHT_ISOSCELES=SQUARE, + PARALLELOGRAMMIC=OBLIQUE, + CENTERED_RHOMBIC=RECTANGULAR, + RIGHT_TRIANGULAR=RECTANGULAR, + CENTERED_RECTANGULAR=RHOMBIC, + ISOSCELE_TRIANGULAR=RHOMBIC, + RIGHT_ISOSCELE_TRIANGULAR=SQUARE, + HEXAGONAL=EQUILATERAL_TRIANGULAR +} LatticeType; + + +/* + * Lagrange-Gauss reduction of a 2D basis. + * The output shall satisfy |out1| <= |out2| <= |out2 - out1| + */ +void l2d_reduceBasis(cart2_t in1, cart2_t in2, cart2_t *out1, cart2_t *out2); + +/* + * This gives the "ordered shortest triple" of base vectors (each pair from the triple + * is a base) and there may not be obtuse angle between o1, o2 and between o2, o3 + */ +void l2d_shortestBase3(cart2_t i1, cart2_t i2, cart2_t *o1, cart2_t *o2, cart2_t *o3); + +/* + * TODO doc + * return value is 4 or 6. + */ +int l2d_shortestBase46(cart2_t i1, cart2_t i2, cart2_t *o1, cart2_t *o2, cart2_t *o3, cart2_t *o4, cart2_t *o5, cart2_t *o6, double rtol); +// variant +int l2d_shortestBase46_arr(cart2_t i1, cart2_t i2, cart2_t *oarr, double rtol); + +// Determines whether angle between inputs is obtuse +bool l2d_is_obtuse(cart2_t i1, cart2_t i2, double rtol); + +/* + * Given two basis vectors, returns 2D Bravais lattice type. + */ +LatticeType l2d_classifyLattice(cart2_t b1, cart2_t b2, double rtol); + +// Other functions in lattices2d.py: TODO? +// range2D() +// generateLattice() +// generateLatticeDisk() +// cutWS() +// filledWS() +// change_basis() + +/* + * Given basis vectors, returns the corners of the Wigner-Seits unit cell (W1, W2, -W1, W2) + * for rectangular and square lattice or (w1, w2, w3, -w1, -w2, -w3) otherwise. + */ +int l2d_cellCornersWS(cart2_t i1, cart2_t i2, cart2_t *o1, cart2_t *o2, cart2_t *o3, cart2_t *o4, cart2_t *o5, cart2_t *o6, double rtol); +// variant +int l2d_cellCornersWS_arr(cart2_t i1, cart2_t i2, cart2_t *oarr, double rtol); + +// Reciprocal bases +void l2d_reciprocalBasis1(cart2_t b1, cart2_t b2, cart2_t *rb1, cart2_t *rb2); +void l2D_reciprocalBasis2pi(cart2_t b1, cart2_t b2, cart2_t *rb1, cart2_t *rb2); + + +/* + * THE MORE OR LESS OK PART + * ======================== + */ + /* * General set of points ordered by the r-coordinate. * Typically, this will include all lattice inside a certain circle. @@ -57,6 +150,21 @@ static inline double points2d_rordered_get_r(const points2d_rordered_t *ps, int return ps->rs[r_order]; } + +ptrdiff_t points2d_rordered_locate_r(const points2d_rordered_t *, double r); + +// returns a "view" (does not copy any of the arrays) +// -- DO NOT FREE orig BEFORE THE END OF SCOPE OF THE RESULT +points2d_rordered_t points2d_rordered_annulus(const points2d_rordered_t *orig, double minr, bool minr_inc, + double maxr, bool maxr_inc); + + + +/* + * THE UGLY PART + * ============= + */ + /* * EQUILATERAL TRIANGULAR LATTICE */ diff --git a/qpms/lattices2d.c b/qpms/lattices2d.c index 7a9cda7..f161a48 100644 --- a/qpms/lattices2d.c +++ b/qpms/lattices2d.c @@ -42,7 +42,43 @@ points2d_rordered_t *points2d_rordered_scale(const points2d_rordered_t *orig, co return p; } +ptrdiff_t points2d_rordered_locate_r(const points2d_rordered_t *p, const double r) { + //if(p->r_rs[0] > r) + // return -1; + //if(p->r_rs[p->nrs-1] < r) + // return p->nrs; + ptrdiff_t lo = 0, hi = p->nrs-1, piv; + while(lo < hi) { + piv = (lo + hi + 1) / 2; + if(p->rs[piv] > r) // the result will be less or equal + hi = piv - 1; + else + lo = piv; + } + return lo; +} +points2d_rordered_t points2d_rordered_annulus(const points2d_rordered_t *orig, + double minr, bool inc_minr, double maxr, bool inc_maxr) { + points2d_rordered_t p; + ptrdiff_t imin, imax; + imin = points2d_rordered_locate_r(orig, minr); + imax = points2d_rordered_locate_r(orig, maxr); + if (!inc_minr && (orig->rs[imin] <= minr)) ++imin; + if (!inc_maxr && (orig->rs[imax] >= maxr)) --imax; + if (imax < imin) { // it's empty + p.nrs = 0; + p.base = NULL; + p.rs = NULL; + p.r_offsets = NULL; + } else { + p.base = orig->base; + p.nrs = imax - imin + 1; + p.rs = orig->rs + imin; + p.r_offsets = orig->r_offsets + imin; + } + return p; +} /* * EQUILATERAL TRIANGULAR LATTICE diff --git a/tests/lattice/2d_triangular_lattice.c b/tests/lattice/2d_triangular_lattice.c index 1d36f7b..7d9ef74 100644 --- a/tests/lattice/2d_triangular_lattice.c +++ b/tests/lattice/2d_triangular_lattice.c @@ -32,6 +32,16 @@ int main() { dump_points2d_rordered(&(g->ps), "triang_h_s160_20.out"); p = points2d_rordered_scale(&(g->ps), -.2); dump_points2d_rordered(p, "triang_h_s160_20_scaled.out"); + { + points2d_rordered_t v = points2d_rordered_annulus(p, 15, true, 17, true); + dump_points2d_rordered(&v, "triang_h_ann_i15-i17.out"); + v = points2d_rordered_annulus(p, 15, false, 17, true); + dump_points2d_rordered(&v, "triang_h_ann_x15-i17.out"); + v = points2d_rordered_annulus(p, 15, false, 17, false); + dump_points2d_rordered(&v, "triang_h_ann_x15-x17.out"); + v = points2d_rordered_annulus(p, 15, true, 17, false); + dump_points2d_rordered(&v, "triang_h_ann_i15-x17.out"); + } points2d_rordered_free(p); triangular_lattice_gen_free(g); @@ -48,7 +58,7 @@ int main() { dump_points2d_rordered(&(g->ps), "triang_v_minus_s7.out"); p = points2d_rordered_scale(&(g->ps), 1/7.); triangular_lattice_gen_free(g); - dump_points2d_rordered(p, "triang_v_minus_s7_scaled_out"); + dump_points2d_rordered(p, "triang_v_minus_s7_scaled.out"); points2d_rordered_free(p); honeycomb_lattice_gen_t *h = honeycomb_lattice_gen_init_h(1, TRIANGULAR_HORIZONTAL);