lattices.h

Go to the documentation of this file.

 
#ifndef LATTICES_H
#define LATTICES_H
 
#include <math.h>
#include <stdbool.h>
#include <assert.h>
#include <stddef.h>
#include <stdlib.h>
#ifndef M_SQRT3
#define M_SQRT3 1.7320508075688772935274463415058724
#endif
#ifndef M_SQRT3_2
#define M_SQRT3_2 (M_SQRT3/2)
#endif
#ifndef M_1_SQRT3
#define M_1_SQRT3 0.57735026918962576450914878050195746
#endif
 
 
 
/* 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
 *
 */
 
typedef enum LatticeDimensionality {
    LAT1D = 1,
    LAT2D = 2,
    LAT3D = 4,
    SPACE1D = 8,
    SPACE2D = 16,
    SPACE3D = 32,
    LAT_1D_IN_3D = 33,
    LAT_2D_IN_3D = 34,
    LAT_3D_IN_3D = 40,
    // special coordinate arrangements (indicating possible optimisations)
    LAT_ZONLY = 64,
    LAT_XYONLY = 128,
    LAT_1D_IN_3D_ZONLY = 97, // LAT1D | SPACE3D | 64
    LAT_2D_IN_3D_XYONLY = 162 // LAT2D | SPACE3D | 128
} LatticeDimensionality;
 
inline static bool LatticeDimensionality_checkflags(
        LatticeDimensionality a, LatticeDimensionality flags_a_has_to_contain) {
    return ((a & flags_a_has_to_contain) == flags_a_has_to_contain);
}
 
// fuck, I already had had suitable type
#include "vectors.h"
typedef cart2_t point2d;
 
static inline point2d point2d_fromxy(const double x, const double y) {
    point2d p;
    p.x = x;
    p.y = y;
    return p;
}
 
 
int qpms_reduce_lattice_basis(double *b, 
        const size_t bsize, 
        const size_t ndim, 
 
        double delta
        );
 
 
typedef struct PGen {
    const struct PGenClassInfo * /*const*/ c; 
    void *stateData;
} PGen;
 
typedef enum PGenPointFlags {
    PGEN_NOTDONE = 2, 
    PGEN_OLD_R = 1, 
    PGEN_SINGLE_R = 16,
    PGEN_AT_Z = 4, 
    PGEN_AT_XY = 8, 
    PGEN_METHOD_UNAVAILABLE = 2048, 
    PGEN_DONE = 0, 
    PGEN_COORDS_CART1 = QPMS_COORDS_CART1,
    PGEN_COORDS_CART2 = QPMS_COORDS_CART2,
    PGEN_COORDS_CART3 = QPMS_COORDS_CART3,
    PGEN_COORDS_POL = QPMS_COORDS_POL,
    PGEN_COORDS_SPH = QPMS_COORDS_SPH,
    PGEN_COORDS_BITRANGE = QPMS_COORDS_BITRANGE,
} PGenPointFlags;
 
 
typedef struct PGenReturnDataBulk {
    PGenPointFlags flags;
    size_t generated; 
} PGenReturnDataBulk;
 
typedef struct PGenReturnData { 
    PGenPointFlags flags; 
    anycoord_point_t point; 
} PGenReturnData;
 
typedef struct PGenZReturnData {
  PGenPointFlags flags; 
  double point_z; 
} PGenZReturnData;
 
typedef struct PGenPolReturnData {
  PGenPointFlags flags; 
  pol_t point_pol; 
} PGenPolReturnData;
 
typedef struct PGenSphReturnData {
  PGenPointFlags flags; 
  sph_t point_sph; 
} PGenSphReturnData;
 
typedef struct PGenCart2ReturnData {
  PGenPointFlags flags; 
  cart2_t point_cart2; 
} PGenCart2ReturnData;
 
typedef struct PGenCart3ReturnData {
  PGenPointFlags flags; 
  cart3_t point_cart3; 
} PGenCart3ReturnData;
 
// convenience constants for use in the extractor implementations
static const PGenZReturnData PGenZDoneVal = {PGEN_DONE, 0}; 
static const PGenPolReturnData PGenPolDoneVal = {PGEN_DONE, {0,0}}; 
static const PGenSphReturnData PGenSphDoneVal = {PGEN_DONE, {0,0,0}}; 
static const PGenCart2ReturnData PGenCart2DoneVal = {PGEN_DONE, {0,0}};
static const PGenCart3ReturnData PGenCart3DoneVal = {PGEN_DONE, {0,0,0}};
 
 
 
typedef struct PGenClassInfo { // static PGenSph info
    char * const name; // mainly for debugging purposes
    int dimensionality; // lower-dimensional can be converted to higher-D, not vice versa; bit redundant with the following, whatever.
    PGenPointFlags native_point_flags;
    PGenReturnData (*next)(struct PGen *);
    PGenZReturnData (*next_z)(struct PGen *);
    PGenPolReturnData (*next_pol)(struct PGen *); 
    PGenSphReturnData (*next_sph)(struct PGen *);
    PGenCart2ReturnData (*next_cart2)(struct PGen *);
    PGenCart3ReturnData (*next_cart3)(struct PGen *);
    PGenReturnDataBulk (*fetch)(struct PGen *, size_t, anycoord_point_t *);
    PGenReturnDataBulk (*fetch_z)(struct PGen *, size_t, double *);
    PGenReturnDataBulk (*fetch_pol)(struct PGen *, size_t, pol_t *); 
    PGenReturnDataBulk (*fetch_sph)(struct PGen *, size_t, sph_t *);
    PGenReturnDataBulk (*fetch_cart2)(struct PGen *, size_t, cart2_t *);
    PGenReturnDataBulk (*fetch_cart3)(struct PGen *, size_t, cart3_t *);
 
    void (*destructor)(struct PGen *); 
} PGenClassInfo;
 
static inline PGenReturnData PGen_next_nf(struct PGen *g) {
    if (g->c->next)
        return g->c->next(g);
    else {
        PGenReturnData r;
        if (g->c->next_z) {
            PGenZReturnData res = g->c->next_z(g);
            r.flags = res.flags;
            r.point.z = res.point_z;
        } else if (g->c->next_pol) {
            PGenPolReturnData res = g->c->next_pol(g);
            r.flags = res.flags;
            r.point.pol = res.point_pol;
        } else if (g->c->next_cart2) {
            PGenCart2ReturnData res = g->c->next_cart2(g);
            r.flags = res.flags;
            r.point.cart2 = res.point_cart2;
        } else if (g->c->next_sph) {
                PGenSphReturnData res = g->c->next_sph(g);
            r.flags = res.flags;
            r.point.sph = res.point_sph;
        } else if (g->c->next_cart3) {
            PGenCart3ReturnData res = g->c->next_cart3(g);
            r.flags = res.flags;
            r.point.cart3 = res.point_cart3;
        } else 
            r.flags = PGEN_METHOD_UNAVAILABLE;
        return r;
    }
}
 
// Ultimate ugliness
static inline PGenReturnDataBulk PGen_fetch_any(struct PGen *g, size_t nmemb,
            anycoord_point_t *target) {
    if (g->c->fetch)
        return g->c->fetch(g, nmemb, target);
    else if (g->c->fetch_cart3) {
        cart3_t *t2 = (cart3_t*) ((char *) target 
            + nmemb * (sizeof(anycoord_point_t)-sizeof(cart3_t)));
        PGenReturnDataBulk res = g->c->fetch_cart3(g, nmemb, t2);
        for (size_t i = 0; i < nmemb; ++i) 
            target[i].cart3 = t2[i];
        return res;
    } else if (g->c->fetch_sph) {
        sph_t *t2 = (sph_t*) ((char *) target 
            + nmemb * (sizeof(anycoord_point_t)-sizeof(sph_t)));
        PGenReturnDataBulk res = g->c->fetch_sph(g, nmemb, t2);
        for (size_t i = 0; i < nmemb; ++i) 
            target[i].sph = t2[i];
        return res;
    } else if (g->c->fetch_cart2) {
        cart2_t *t2 = (cart2_t*) ((char *) target 
            + nmemb * (sizeof(anycoord_point_t)-sizeof(cart2_t)));
        PGenReturnDataBulk res = g->c->fetch_cart2(g, nmemb, t2);
        for (size_t i = 0; i < nmemb; ++i) 
            target[i].cart2 = t2[i];
        return res;
    } else if (g->c->fetch_pol) {
        pol_t *t2 = (pol_t*) ((char *) target 
            + nmemb * (sizeof(anycoord_point_t)-sizeof(pol_t)));
        PGenReturnDataBulk res = g->c->fetch_pol(g, nmemb, t2);
        for (size_t i = 0; i < nmemb; ++i) 
            target[i].pol = t2[i];
        return res;
    } else if (g->c->fetch_z) {
        double *t2 = (double*) ((char *) target 
            + nmemb * (sizeof(anycoord_point_t)-sizeof(double)));
        PGenReturnDataBulk res = g->c->fetch_z(g, nmemb, t2);
        for (size_t i = 0; i < nmemb; ++i) 
            target[i].z = t2[i];
        return res;
    } else {
        // This is ridiculously inefficient
        PGenReturnDataBulk res = {PGEN_NOTDONE, 0};
        for (res.generated = 0; res.generated < nmemb;
                    ++res.generated) {
            PGenReturnData res1 = PGen_next_nf(g);
            QPMS_ENSURE(!(res1.flags & PGEN_METHOD_UNAVAILABLE),
                "No method found to generate points. The PGenClassInfo"
                " %s is apparently broken.", g->c->name);
            if (res1.flags & PGEN_NOTDONE) {
                target[res.generated] = res1.point;
                // The coordinate system generated by next() must be consistent:
                assert(!res.generated || ((res1.flags & PGEN_COORDS_BITRANGE) == (res.flags & PGEN_COORDS_BITRANGE)));
                res.flags |= res1.flags & PGEN_COORDS_BITRANGE;
            } else {
                res.flags &= ~PGEN_NOTDONE;
                break;
            }
        }
        // Do not guarantee anything for; low priority TODO
        res.flags &= ~(PGEN_OLD_R & PGEN_SINGLE_R);
        return res;
    }
}
 
static inline PGenReturnData PGen_next(struct PGen *g) {
    PGenReturnData res = PGen_next_nf(g);
    if (!(res.flags & PGEN_METHOD_UNAVAILABLE))
        return res;
    else { // Slow if implementation is stupid, but short!
        PGenReturnDataBulk resb = PGen_fetch_any(g, 1, &res.point);
        if (resb.generated) 
            // the | PGEN_NOTDONE may not be needed, but my brain melted
            res.flags = resb.flags | PGEN_NOTDONE;
        else
            res.flags = PGEN_DONE;
        return res;
    }
}
 
static inline PGenReturnDataBulk PGen_fetch_sph(struct PGen *g, 
        size_t nmemb, sph_t *target) {
    if (g->c->fetch_sph)
        return g->c->fetch_sph(g, nmemb, target);
    else {
        anycoord_point_t *tmp;
        QPMS_CRASHING_MALLOC(tmp, sizeof(anycoord_point_t) * nmemb);
        PGenReturnDataBulk res = PGen_fetch_any(g, nmemb, tmp);
        anycoord_arr2something(target, QPMS_COORDS_SPH,
                tmp, res.flags, res.generated);
        free(tmp);
        res.flags = (res.flags & ~QPMS_COORDS_BITRANGE)
            | QPMS_COORDS_SPH;
        return res;
    }
}
 
static inline PGenReturnDataBulk PGen_fetch_cart3(struct PGen *g, 
        size_t nmemb, cart3_t *target) {
    if (g->c->fetch_cart3)
        return g->c->fetch_cart3(g, nmemb, target);
    else {
        anycoord_point_t *tmp;
        QPMS_CRASHING_MALLOC(tmp, sizeof(anycoord_point_t) * nmemb);
        PGenReturnDataBulk res = PGen_fetch_any(g, nmemb, tmp);
        anycoord_arr2something(target, QPMS_COORDS_CART3,
                tmp, res.flags, res.generated);
        free(tmp);
        res.flags = (res.flags & ~QPMS_COORDS_BITRANGE)
            | QPMS_COORDS_CART3;
        return res;
    }
}
 
static inline PGenReturnDataBulk PGen_fetch_pol(struct PGen *g, 
        size_t nmemb, pol_t *target) {
    if (g->c->fetch_pol)
        return g->c->fetch_pol(g, nmemb, target);
    else {
        anycoord_point_t *tmp;
        QPMS_CRASHING_MALLOC(tmp, sizeof(anycoord_point_t) * nmemb);
        PGenReturnDataBulk res = PGen_fetch_any(g, nmemb, tmp);
        anycoord_arr2something(target, QPMS_COORDS_POL,
                tmp, res.flags, res.generated);
        free(tmp);
        res.flags = (res.flags & ~QPMS_COORDS_BITRANGE)
            | QPMS_COORDS_POL;
        return res;
    }
}
 
static inline PGenReturnDataBulk PGen_fetch_cart2(struct PGen *g, 
        size_t nmemb, cart2_t *target) {
    if (g->c->fetch_cart2)
        return g->c->fetch_cart2(g, nmemb, target);
    else {
        anycoord_point_t *tmp;
        QPMS_CRASHING_MALLOC(tmp, sizeof(anycoord_point_t) * nmemb);
        PGenReturnDataBulk res = PGen_fetch_any(g, nmemb, tmp);
        anycoord_arr2something(target, QPMS_COORDS_CART2,
                tmp, res.flags, res.generated);
        free(tmp);
        res.flags = (res.flags & ~QPMS_COORDS_BITRANGE)
            | QPMS_COORDS_CART2;
        return res;
    }
}
 
static inline PGenReturnDataBulk PGen_fetch_z(struct PGen *g, 
        size_t nmemb, double *target) {
    if (g->c->fetch_z)
        return g->c->fetch_z(g, nmemb, target);
    else {
        anycoord_point_t *tmp;
        QPMS_CRASHING_MALLOC(tmp, sizeof(anycoord_point_t) * nmemb);
        PGenReturnDataBulk res = PGen_fetch_any(g, nmemb, tmp);
        anycoord_arr2something(target, QPMS_COORDS_CART1,
                tmp, res.flags, res.generated);
        free(tmp);
        res.flags = (res.flags & ~QPMS_COORDS_BITRANGE)
            | QPMS_COORDS_CART1;
        return res;
    }
}
 
static inline void PGen_destroy(PGen *g) {
    g->c->destructor(g);
    assert(g->stateData == NULL); // this should be done by the destructor
}
 
static inline PGenZReturnData PGen_next_z(PGen *g) {
    if (g->c->next_z)
        return g->c->next_z(g);
    else { // Super-slow generic fallback.
        PGenReturnData res = PGen_next(g);
        if (res.flags & PGEN_NOTDONE) {
            PGenZReturnData r;
            r.point_z = anycoord2cart1(res.point, res.flags);
            r.flags = (res.flags & ~QPMS_COORDS_BITRANGE)
                | QPMS_COORDS_CART1;
            return r;
        } else 
            return PGenZDoneVal;
    }
}
 
static inline PGenSphReturnData PGen_next_sph(PGen *g) {
    if (g->c->next_sph) 
        return g->c->next_sph(g);
    else { // Super-slow generic fallback.
        PGenReturnData res = PGen_next(g);
        if (res.flags & PGEN_NOTDONE) {
            PGenSphReturnData r;
            r.point_sph = anycoord2sph(res.point, res.flags);
            r.flags = (res.flags & ~QPMS_COORDS_BITRANGE)
                | QPMS_COORDS_SPH;
            return r;
        } else 
            return PGenSphDoneVal;
    }
}
 
static inline PGenPolReturnData PGen_next_pol(PGen *g) {
    if (g->c->next_pol) 
        return g->c->next_pol(g);
    else { // Super-slow generic fallback.
        PGenReturnData res = PGen_next(g);
        if (res.flags & PGEN_NOTDONE) {
            PGenPolReturnData r;
            r.point_pol = anycoord2pol(res.point, res.flags);
            r.flags = (res.flags & ~QPMS_COORDS_BITRANGE)
                | QPMS_COORDS_POL;
            return r;
        } else 
            return PGenPolDoneVal;
    }
}
 
static inline PGenCart3ReturnData PGen_next_cart3(PGen *g) {
    if (g->c->next_cart3) 
        return g->c->next_cart3(g);
    else { // Super-slow generic fallback.
        PGenReturnData res = PGen_next(g);
        if (res.flags & PGEN_NOTDONE) {
            PGenCart3ReturnData r;
            r.point_cart3 = anycoord2cart3(res.point, res.flags);
            r.flags = (res.flags & ~QPMS_COORDS_BITRANGE)
                | QPMS_COORDS_CART3;
            return r;
        } else 
            return PGenCart3DoneVal;
    }
}
 
static inline PGenCart2ReturnData PGen_next_cart2(PGen *g) {
    if (g->c->next_cart2) 
        return g->c->next_cart2(g);
    else { // Super-slow generic fallback.
        PGenReturnData res = PGen_next(g);
        if (res.flags & PGEN_NOTDONE) {
            PGenCart2ReturnData r;
            r.point_cart2 = anycoord2cart2(res.point, res.flags);
            r.flags = (res.flags & ~QPMS_COORDS_BITRANGE)
                | QPMS_COORDS_CART2;
            return r;
        } else 
            return PGenCart2DoneVal;
    }
}
 
#if 0
static inline PGenReturnDataBulk PGen_fetch(PGen *g, size_t n, anycoord_point_t *arr){
    if (g->c->fetch)
        return g->c->fetch(g, n, arr);
    else abort();
}
#endif
 
static inline bool PGenSph_notDone(PGenSphReturnData data) {
    return data.flags & PGEN_NOTDONE ? true : false;
}
static inline bool PGenCart3_notDone(PGenCart3ReturnData data) {
    return data.flags & PGEN_NOTDONE ? true : false;
}
 
static inline PGenCart3ReturnData PGenReturnDataConv_sph_cart3(PGenSphReturnData sphdata){
    PGenCart3ReturnData c3data;
    c3data.flags = sphdata.flags;
    c3data.point_cart3 = sph2cart(sphdata.point_sph);
    return c3data; 
}
 
static inline PGenSphReturnData PGenReturnDataConv_cart3_sph(PGenCart3ReturnData c){
    PGenSphReturnData s;
    s.flags = c.flags;
    s.point_sph = cart2sph(c.point_cart3);
    return s; 
}
 
/*
 * Some basic lattice generators implementing the abstract interface above (implemented in latticegens.c).
 */
 
extern const PGenClassInfo PGen_FromPoint2DArray; // TODO Do I even need this to be declared here?
PGen PGen_FromPoints2DArray_new(const point2d *points, size_t len);
 
extern const PGenClassInfo PGen_1D;
typedef enum PGen_1D_incrementDirection{
    //PGEN_1D_POSITIVE_INC, // not implemented
    //PGEN_1D_NEGATIVE_INC, // not implemented
    PGEN_1D_INC_FROM_ORIGIN,
    PGEN_1D_INC_TOWARDS_ORIGIN
} PGen_1D_incrementDirection;
PGen PGen_1D_new_minMaxR(double period, 
    double offset, 
    double minR, 
    bool inc_minR, 
    double maxR, 
    bool inc_maxR, 
    PGen_1D_incrementDirection incdir 
);
 
extern const PGenClassInfo PGen_xyWeb;
PGen PGen_xyWeb_new(cart2_t b1, cart2_t b2, double rtol, cart2_t offset, 
        double minR, bool inc_minR, double maxR, bool inc_maxR);
 
size_t PGen_xyWeb_sizecap(cart2_t b1, cart2_t b2, double rtol, cart2_t offset,
        double minR, bool inc_minR, double maxR, bool inc_maxR);
 
 
extern const PGenClassInfo PGen_LatticeRadialHeap2D;
extern const PGenClassInfo PGen_LatticeRadialHeap3D;
PGen PGen_LatticeRadialHeap2D_new(cart2_t b1, cart2_t b2, cart2_t offset, 
        double minR, bool inc_minR, double maxR, bool inc_maxR);
PGen PGen_LatticeRadialHeap3D_new(const cart3_t *b1, const cart3_t *b2, const cart3_t *b3,
            const cart3_t *offset, double minR, bool inc_minR, double maxR, bool inc_maxR);
 
 
extern const PGenClassInfo PGen_shifted;
PGen PGen_shifted_new(PGen orig, cart3_t shift);
 
/*
 * 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
 
// Bravais lattice types
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
} LatticeType2;
 
#if 0
// Wallpaper groups
typedef enum {
    TODO
} SpaceGroup2;
#endif
 
// Just for detecting the right angles (needed for generators).
typedef enum {
    NOT_ORTHOGONAL = 0,
    ORTHOGONAL_01 = 1,
    ORTHOGONAL_12 = 2,
    ORTHOGONAL_02 = 4
} LatticeFlags;
 
 
 
/*
 * 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 one uses LLL reduction.
void l3d_reduceBasis(const cart3_t in[3], cart3_t out[3]);
 
/* 
 * 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_r(cart2_t i1, cart2_t i2, double rtol);
bool l2d_is_obtuse(cart2_t i1, cart2_t i2);
 
/* 
 * Given two basis vectors, returns 2D Bravais lattice type.
 */
LatticeType2 l2d_classifyLattice(cart2_t b1, cart2_t b2, double rtol);
 
// Detects right angles.
LatticeFlags l2d_detectRightAngles(cart2_t b1, cart2_t b2, double rtol);
LatticeFlags l3d_detectRightAngles(const cart3_t basis[3], 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; returns 0 on success, possibly a non-zero if b1 and b2 are parallel
int l2d_reciprocalBasis1(cart2_t b1, cart2_t b2, cart2_t *rb1, cart2_t *rb2);
int l2d_reciprocalBasis2pi(cart2_t b1, cart2_t b2, cart2_t *rb1, cart2_t *rb2);
// 3D reciprocal bases; returns (direct) unit cell volume with possible sign. Assumes direct lattice basis already reduced.
double l3d_reciprocalBasis1(const cart3_t direct_basis[3], cart3_t reciprocal_basis[3]);
double l3d_reciprocalBasis2pi(const cart3_t direct_basis[3], cart3_t reciprocal_basis[3]);
 
double l2d_unitcell_area(cart2_t b1, cart2_t b2);
 
// returns the radius of inscribed circle of a hexagon (or rectangle/square if applicable) created by the shortest base triple
double l2d_hexWebInCircleRadius(cart2_t b1, cart2_t b2);
 
/*
 * 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.
 * This structure is internally used by the "lattice generators" below.
 * It does not have its memory management of its own, as it is handled
 * by the "generators". For everything except the generators,
 * this structure shall be read-only.
 */
typedef struct {
    size_t nrs; // number of different radii
    double *rs; // the radii; of length nrs (largest contained radius == rs[nrs-1])
    point2d *base; 
    ptrdiff_t *r_offsets; // of length nrs+1 (using relative offsets due to possible realloc's)
    // the jth point of i-th radius is base[r_offsets[i]+j] or using the inline below..
    /* // redundand (therefore removed) members
         * point2d *points; // redundant as it is the same as points_at_r[0]
     * size_t npoints;  // redundant as it is the same as points_at_r[nrs]-points_at_r[0]
     */ 
} points2d_rordered_t;
 
 
// sorts arbitrary points and creates points2d_rordered_t
points2d_rordered_t *points2d_rordered_frompoints(const point2d *orig_base,
            size_t nmemb, double rtol, double atol);
 
// returns a copy but shifted by a constant (actually in a stupid way, but whatever)
points2d_rordered_t *points2d_rordered_shift(const points2d_rordered_t *orig,
        point2d shift, double rtol, double atol);
 
// returns a copy but scaled by a factor
points2d_rordered_t *points2d_rordered_scale(const points2d_rordered_t *orig,
        double factor);
 
 
/* The destructor: use only for results of
 *  - points2D_rordered_frompoints,
 *  - points2d_rordered_shift,
 *  - points2d_rordered_scale.
 */
void points2d_rordered_free(points2d_rordered_t *);
 
static inline point2d points2d_rordered_get_point(const points2d_rordered_t *ps, int r_order, int i) {
    assert(i >= 0);
    assert(r_order < ps->nrs);
    assert(i < (ps->r_offsets[r_order+1] - ps->r_offsets[r_order]));
    return ps->base[ps->r_offsets[r_order] + i];
}
 
static inline double points2d_rordered_get_r(const points2d_rordered_t *ps, int r_order) {
    assert(r_order < ps->nrs);
    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);
 
 
#if 0
double qpms_emptylattice2_mode_nth(
                cart2_t b1_rec, 
                cart2_t b2_rec, 
                double rtol, 
                cart2_t k, 
                double wave_speed, 
                size_t N 
                );
 
void qpms_emptylattice2_modes_maxindex(
                double target_freqs[], 
                cart2_t b1_rec, 
                cart2_t b2_rec, 
                double rtol, 
                cart2_t k, 
                double wave_speed, 
                size_t maxindex 
                );
#endif
 
 
size_t qpms_emptylattice2_modes_maxfreq(
                double **target_freqs,
                cart2_t b1_rec, 
                cart2_t b2_rec, 
                double rtol, 
                cart2_t k, 
                double wave_speed, 
                double maxfreq 
                );
 
void qpms_emptylattice2_modes_nearest(
                double target[2], 
                cart2_t b1_rec, 
                cart2_t b2_rec, 
                double rtol, 
                cart2_t k, 
                double wave_speed, 
                double omega 
                );
 
 
 
/*
 * THE UGLY PART
 * =============
 */
 
/* 
 * EQUILATERAL TRIANGULAR LATTICE
 */
 
typedef enum {
    TRIANGULAR_VERTICAL, // there is a lattice base vector parallel to the y-axis
    TRIANGULAR_HORIZONTAL // there is a lattice base vector parallel to the x-axis
} TriangularLatticeOrientation;
 
static inline TriangularLatticeOrientation reverseTriangularLatticeOrientation(TriangularLatticeOrientation o){
    switch(o) {
        case TRIANGULAR_VERTICAL:
            return TRIANGULAR_HORIZONTAL;
            break;
        case TRIANGULAR_HORIZONTAL:
            return TRIANGULAR_VERTICAL;
            break;
        default:
            abort();
    }
    abort();
}
 
// implementation data structures; not needed in the header file
typedef struct triangular_lattice_gen_privstuff_t triangular_lattice_gen_privstuff_t;
 
typedef struct {
    // public:
    points2d_rordered_t ps;
    TriangularLatticeOrientation orientation;
    double a; // lattice vector length
    
    // not sure if needed:
    bool includes_origin;
 
    // Denotes an offset of the "origin" point; meaning step hexshift * a / sqrt(2) upwards
    // or leftwards for the horizontal or vertical orientations, respectively.
    int hexshift; 
 
    // private:
    triangular_lattice_gen_privstuff_t *priv;
 
} triangular_lattice_gen_t;
 
triangular_lattice_gen_t *triangular_lattice_gen_init(double a, TriangularLatticeOrientation ori, bool include_origin,
        int halfoffset); 
const points2d_rordered_t * triangular_lattice_gen_getpoints(const triangular_lattice_gen_t *g);
int triangular_lattice_gen_extend_to_r(triangular_lattice_gen_t *g, double r);
int triangular_lattice_gen_extend_to_steps(triangular_lattice_gen_t *g, int maxsteps);
void triangular_lattice_gen_free(triangular_lattice_gen_t *g);
 
 
/*
 * HONEYCOMB LATTICE
 */
 
typedef struct {
    // public:
    points2d_rordered_t ps;
    TriangularLatticeOrientation orientation;
    double a;
    double h;
 
    // private:
    triangular_lattice_gen_t *tg;
} honeycomb_lattice_gen_t;
 
honeycomb_lattice_gen_t *honeycomb_lattice_gen_init_h(double h, TriangularLatticeOrientation ori);
honeycomb_lattice_gen_t *honeycomb_lattice_gen_init_a(double a, TriangularLatticeOrientation ori);
int honeycomb_lattice_gen_extend_to_steps(honeycomb_lattice_gen_t *g, int maxsteps);
int honeycomb_lattice_gen_extend_to_r(honeycomb_lattice_gen_t *g, double r);
void honeycomb_lattice_gen_free(honeycomb_lattice_gen_t *g);
 
#endif // LATTICES_H