diff --git a/qpms/latticegens.c b/qpms/latticegens.c index 9bc2ac0..f556b4c 100644 --- a/qpms/latticegens.c +++ b/qpms/latticegens.c @@ -432,8 +432,12 @@ PGenCart2ReturnData PGen_xyWeb_next_cart2(PGen *g) { if (s->layer <= s->last_layer) { const cart2_t thePoint = cart2_add(s->offset, cart2_add(cart2_scale(s->i, s->b1), cart2_scale(s->j, s->b2))); - if(s->layer == 0) // origin is unique, proceed with next layer + if(s->layer == 0) { // origin is unique, proceed with next layer ++s->layer; + s->phase = 0; + s->i = s->layer; + s->j = 0; + } else if(s->lt & (SQUARE | RECTANGULAR)) { // rectangular or square lattice, four perimeters switch(s->phase) { diff --git a/qpms/lattices2d.c b/qpms/lattices2d.c index 6f4e1b4..8bf6778 100644 --- a/qpms/lattices2d.c +++ b/qpms/lattices2d.c @@ -733,6 +733,39 @@ int l2d_shortestBase46(const cart2_t i1, const cart2_t i2, cart2_t *o1, cart2_t } } +/* + * Given two basis vectors, returns 2D Bravais lattice type. + */ +LatticeType2 l2d_classifyLattice(cart2_t b1, cart2_t b2, double rtol) +{ + l2d_reduceBasis(b1, b2, &b1, &b2); + cart2_t b3 = cart2_substract(b2, b1); + double b1s = cart2_normsq(b1), b2s = cart2_normsq(b2), b3s = cart2_normsq(b3); + double eps = rtol * (b2s + b1s); //FIXME what should eps be? + // avoid obtuse angle between b1 and b2. TODO this should be yet tested + // TODO use is_obtuse here? + if (b3s - b2s - b1s > eps) { + b3 = b2; + b2 = cart2_add(b2, b1); + // N.B. now the assumption |b3| >= |b2| is no longer valid + // b3 = cart2_substract(b2, b1) + b2s = cart2_normsq(b2); + b3s = cart2_normsq(b3); + } + if (fabs(b2s-b1s) < eps || fabs(b2s - b3s) < eps) { // isoscele + if (fabs(b3s-b1s) < eps) + return EQUILATERAL_TRIANGULAR; + else if (fabs(b3s - 2*b1s)) + return SQUARE; + else + return RHOMBIC; + } else if (fabs(b3s-b2s-b1s) < eps) + return RECTANGULAR; + else + return OBLIQUE; +} + + # if 0 // variant int l2d_shortestBase46_arr(cart2_t i1, cart2_t i2, cart2_t *oarr, double rtol); @@ -740,10 +773,6 @@ 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); -/* - * 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() diff --git a/tests/test_latticegenxyweb.c b/tests/test_latticegenxyweb.c new file mode 100644 index 0000000..2c98b3d --- /dev/null +++ b/tests/test_latticegenxyweb.c @@ -0,0 +1,54 @@ +//c99 -o test_latticegenxyweb -ggdb -Wall -I ../ test_latticegenxyweb.c ../qpms/latticegens.c -lm +#define QPMS_VECTORS_NICE_TRANSFORMATIONS +#include +#include + +void fprint_PGenCart2ReturnData(FILE *f, PGenCart2ReturnData d) { + cart2_t c = d.point_cart2; + fprintf(f, "%g\t%g\tflags: %#xd\n", + c.x, c.y, d.flags); +} + +void dump_PGenCart2(char *filename, PGen *g) { + FILE *out = fopen(filename, "w"); + PGenCart2ReturnData d; + while ((d = PGen_next_cart2(g)).flags & PGEN_NOTDONE) { + fprint_PGenCart2ReturnData(out, d); + } +} + +#if 0 +void print_PGenSphReturnData(PGenSphReturnData d) { + sph_t s = d.point_sph; + cart3_t c = sph2cart(s); + printf("sph: (%g, %g π, %g π), cart: (%g, %g, %g), flags: %#xd\n", + s.r, s.theta / M_PI, s.phi / M_PI, c.x, c.y, c.z, d.flags); +} + +void dump_PGenSph(PGen *g) { + PGenSphReturnData d; + do { + d = PGen_next_sph(g); + print_PGenSphReturnData(d); + } while (d.flags & PGEN_NOTDONE); +} +#endif + +#define DO_AND_PRINT(label, s) printf(#label ":\n" #s "\n"); s ; + +int main(int argc, char **argv) { + PGen g; + cart2_t b1 = {1.1,0}, b2 = {0,-1.3}, offset = {0.1, -0.1}; + g = PGen_xyWeb_new(b1, b2, 1e-13, offset, 0, true, 5, false); + dump_PGenCart2("rect1.xydump", &g); + g = PGen_xyWeb_new(b1, b2, 1e-13, offset, 5, true, 8, false); + dump_PGenCart2("rect2.xydump", &g); + b1.x = 1.342, b1.y = 4.3121; b2.x = -1.83, b2.y = 1.4; + g = PGen_xyWeb_new(b1, b2, 1e-13, offset, 0, true, 8, false); + dump_PGenCart2("oblique1.xydump", &g); + + return 0; + + +} +