QPMS
Electromagnetic multiple scattering library and toolkit.
lattices.h
Go to the documentation of this file.
1 
5 #ifndef LATTICES_H
6 #define LATTICES_H
7 
8 #include <math.h>
9 #include <stdbool.h>
10 #include <assert.h>
11 #include <stddef.h>
12 #include <stdlib.h>
13 #ifndef M_SQRT3
14 #define M_SQRT3 1.7320508075688772935274463415058724
15 #endif
16 #ifndef M_SQRT3_2
17 #define M_SQRT3_2 (M_SQRT3/2)
18 #endif
19 #ifndef M_1_SQRT3
20 #define M_1_SQRT3 0.57735026918962576450914878050195746
21 #endif
22 
23 
24 
25 /* IMPORTANT TODO
26  * ==============
27  *
28  * The current content of this part (and the implementation) is extremely ugly.
29  * When I have some time, I have to rewrite this in the style of lattices2d.py
30  *
31  */
32 
33 typedef enum LatticeDimensionality {
34  LAT1D = 1,
35  LAT2D = 2,
36  LAT3D = 4,
37  SPACE1D = 8,
38  SPACE2D = 16,
39  SPACE3D = 32,
40  LAT_1D_IN_3D = 33,
41  LAT_2D_IN_3D = 34,
42  LAT_3D_IN_3D = 40,
43  // special coordinate arrangements (indicating possible optimisations)
44  LAT_ZONLY = 64,
45  LAT_XYONLY = 128,
46  LAT_1D_IN_3D_ZONLY = 97, // LAT1D | SPACE3D | 64
47  LAT_2D_IN_3D_XYONLY = 162 // LAT2D | SPACE3D | 128
48 } LatticeDimensionality;
49 
50 inline static bool LatticeDimensionality_checkflags(
51  LatticeDimensionality a, LatticeDimensionality flags_a_has_to_contain) {
52  return ((a & flags_a_has_to_contain) == flags_a_has_to_contain);
53 }
54 
55 // fuck, I already had had suitable type
56 #include "vectors.h"
57 typedef cart2_t point2d;
58 
59 static inline point2d point2d_fromxy(const double x, const double y) {
60  point2d p;
61  p.x = x;
62  p.y = y;
63  return p;
64 }
65 
67 
73 int qpms_reduce_lattice_basis(double *b,
74  const size_t bsize,
75  const size_t ndim,
77 
79  double delta
80  );
81 
83 
120 typedef struct PGen {
122  const struct PGenClassInfo * /*const*/ c;
124  void *stateData;
126 
127 typedef enum PGenPointFlags {
147  PGEN_AT_Z = 4,
150  PGEN_DONE = 0,
151  PGEN_COORDS_CART1 = QPMS_COORDS_CART1,
152  PGEN_COORDS_CART2 = QPMS_COORDS_CART2,
153  PGEN_COORDS_CART3 = QPMS_COORDS_CART3,
154  PGEN_COORDS_POL = QPMS_COORDS_POL,
155  PGEN_COORDS_SPH = QPMS_COORDS_SPH,
156  PGEN_COORDS_BITRANGE = QPMS_COORDS_BITRANGE,
158 
159 
161 typedef struct PGenReturnDataBulk {
164  size_t generated;
166 
168 typedef struct PGenReturnData {
172 
174 typedef struct PGenZReturnData {
176  double point_z;
178 
180 typedef struct PGenPolReturnData {
184 
186 typedef struct PGenSphReturnData {
190 
192 typedef struct PGenCart2ReturnData {
196 
198 typedef struct PGenCart3ReturnData {
202 
203 // convenience constants for use in the extractor implementations
204 static const PGenZReturnData PGenZDoneVal = {PGEN_DONE, 0};
205 static const PGenPolReturnData PGenPolDoneVal = {PGEN_DONE, {0,0}};
206 static const PGenSphReturnData PGenSphDoneVal = {PGEN_DONE, {0,0,0}};
207 static const PGenCart2ReturnData PGenCart2DoneVal = {PGEN_DONE, {0,0}};
208 static const PGenCart3ReturnData PGenCart3DoneVal = {PGEN_DONE, {0,0,0}};
209 
210 
212 
231 typedef struct PGenClassInfo { // static PGenSph info
232  char * const name; // mainly for debugging purposes
233  int dimensionality; // lower-dimensional can be converted to higher-D, not vice versa; bit redundant with the following, whatever.
237  PGenReturnData (*next)(struct PGen *);
251  PGenReturnDataBulk (*fetch_z)(struct PGen *, size_t, double *);
253  PGenReturnDataBulk (*fetch_pol)(struct PGen *, size_t, pol_t *);
255  PGenReturnDataBulk (*fetch_sph)(struct PGen *, size_t, sph_t *);
257  PGenReturnDataBulk (*fetch_cart2)(struct PGen *, size_t, cart2_t *);
259  PGenReturnDataBulk (*fetch_cart3)(struct PGen *, size_t, cart3_t *);
261 
264  void (*destructor)(struct PGen *);
266 
268 static inline PGenReturnData PGen_next_nf(struct PGen *g) {
269  if (g->c->next)
270  return g->c->next(g);
271  else {
272  PGenReturnData r;
273  if (g->c->next_z) {
274  PGenZReturnData res = g->c->next_z(g);
275  r.flags = res.flags;
276  r.point.z = res.point_z;
277  } else if (g->c->next_pol) {
278  PGenPolReturnData res = g->c->next_pol(g);
279  r.flags = res.flags;
280  r.point.pol = res.point_pol;
281  } else if (g->c->next_cart2) {
282  PGenCart2ReturnData res = g->c->next_cart2(g);
283  r.flags = res.flags;
284  r.point.cart2 = res.point_cart2;
285  } else if (g->c->next_sph) {
286  PGenSphReturnData res = g->c->next_sph(g);
287  r.flags = res.flags;
288  r.point.sph = res.point_sph;
289  } else if (g->c->next_cart3) {
290  PGenCart3ReturnData res = g->c->next_cart3(g);
291  r.flags = res.flags;
292  r.point.cart3 = res.point_cart3;
293  } else
295  return r;
296  }
297 }
298 
300 // Ultimate ugliness
301 static inline PGenReturnDataBulk PGen_fetch_any(struct PGen *g, size_t nmemb,
302  anycoord_point_t *target) {
303  if (g->c->fetch)
304  return g->c->fetch(g, nmemb, target);
305  else if (g->c->fetch_cart3) {
306  cart3_t *t2 = (cart3_t*) ((char *) target
307  + nmemb * (sizeof(anycoord_point_t)-sizeof(cart3_t)));
308  PGenReturnDataBulk res = g->c->fetch_cart3(g, nmemb, t2);
309  for (size_t i = 0; i < nmemb; ++i)
310  target[i].cart3 = t2[i];
311  return res;
312  } else if (g->c->fetch_sph) {
313  sph_t *t2 = (sph_t*) ((char *) target
314  + nmemb * (sizeof(anycoord_point_t)-sizeof(sph_t)));
315  PGenReturnDataBulk res = g->c->fetch_sph(g, nmemb, t2);
316  for (size_t i = 0; i < nmemb; ++i)
317  target[i].sph = t2[i];
318  return res;
319  } else if (g->c->fetch_cart2) {
320  cart2_t *t2 = (cart2_t*) ((char *) target
321  + nmemb * (sizeof(anycoord_point_t)-sizeof(cart2_t)));
322  PGenReturnDataBulk res = g->c->fetch_cart2(g, nmemb, t2);
323  for (size_t i = 0; i < nmemb; ++i)
324  target[i].cart2 = t2[i];
325  return res;
326  } else if (g->c->fetch_pol) {
327  pol_t *t2 = (pol_t*) ((char *) target
328  + nmemb * (sizeof(anycoord_point_t)-sizeof(pol_t)));
329  PGenReturnDataBulk res = g->c->fetch_pol(g, nmemb, t2);
330  for (size_t i = 0; i < nmemb; ++i)
331  target[i].pol = t2[i];
332  return res;
333  } else if (g->c->fetch_z) {
334  double *t2 = (double*) ((char *) target
335  + nmemb * (sizeof(anycoord_point_t)-sizeof(double)));
336  PGenReturnDataBulk res = g->c->fetch_z(g, nmemb, t2);
337  for (size_t i = 0; i < nmemb; ++i)
338  target[i].z = t2[i];
339  return res;
340  } else {
341  // This is ridiculously inefficient
342  PGenReturnDataBulk res = {PGEN_NOTDONE, 0};
343  for (res.generated = 0; res.generated < nmemb;
344  ++res.generated) {
345  PGenReturnData res1 = PGen_next_nf(g);
347  "No method found to generate points. The PGenClassInfo"
348  " %s is apparently broken.", g->c->name);
349  if (res1.flags & PGEN_NOTDONE) {
350  target[res.generated] = res1.point;
351  // The coordinate system generated by next() must be consistent:
352  assert(!res.generated || ((res1.flags & PGEN_COORDS_BITRANGE) == (res.flags & PGEN_COORDS_BITRANGE)));
353  res.flags |= res1.flags & PGEN_COORDS_BITRANGE;
354  } else {
355  res.flags &= ~PGEN_NOTDONE;
356  break;
357  }
358  }
359  // Do not guarantee anything for; low priority TODO
360  res.flags &= ~(PGEN_OLD_R & PGEN_SINGLE_R);
361  return res;
362  }
363 }
364 
366 static inline PGenReturnData PGen_next(struct PGen *g) {
367  PGenReturnData res = PGen_next_nf(g);
368  if (!(res.flags & PGEN_METHOD_UNAVAILABLE))
369  return res;
370  else { // Slow if implementation is stupid, but short!
371  PGenReturnDataBulk resb = PGen_fetch_any(g, 1, &res.point);
372  if (resb.generated)
373  // the | PGEN_NOTDONE may not be needed, but my brain melted
374  res.flags = resb.flags | PGEN_NOTDONE;
375  else
376  res.flags = PGEN_DONE;
377  return res;
378  }
379 }
380 
382 static inline PGenReturnDataBulk PGen_fetch_sph(struct PGen *g,
383  size_t nmemb, sph_t *target) {
384  if (g->c->fetch_sph)
385  return g->c->fetch_sph(g, nmemb, target);
386  else {
387  anycoord_point_t *tmp;
388  QPMS_CRASHING_MALLOC(tmp, sizeof(anycoord_point_t) * nmemb);
389  PGenReturnDataBulk res = PGen_fetch_any(g, nmemb, tmp);
391  tmp, res.flags, res.generated);
392  free(tmp);
393  res.flags = (res.flags & ~QPMS_COORDS_BITRANGE)
394  | QPMS_COORDS_SPH;
395  return res;
396  }
397 }
398 
400 static inline PGenReturnDataBulk PGen_fetch_cart3(struct PGen *g,
401  size_t nmemb, cart3_t *target) {
402  if (g->c->fetch_cart3)
403  return g->c->fetch_cart3(g, nmemb, target);
404  else {
405  anycoord_point_t *tmp;
406  QPMS_CRASHING_MALLOC(tmp, sizeof(anycoord_point_t) * nmemb);
407  PGenReturnDataBulk res = PGen_fetch_any(g, nmemb, tmp);
409  tmp, res.flags, res.generated);
410  free(tmp);
411  res.flags = (res.flags & ~QPMS_COORDS_BITRANGE)
413  return res;
414  }
415 }
416 
418 static inline PGenReturnDataBulk PGen_fetch_pol(struct PGen *g,
419  size_t nmemb, pol_t *target) {
420  if (g->c->fetch_pol)
421  return g->c->fetch_pol(g, nmemb, target);
422  else {
423  anycoord_point_t *tmp;
424  QPMS_CRASHING_MALLOC(tmp, sizeof(anycoord_point_t) * nmemb);
425  PGenReturnDataBulk res = PGen_fetch_any(g, nmemb, tmp);
427  tmp, res.flags, res.generated);
428  free(tmp);
429  res.flags = (res.flags & ~QPMS_COORDS_BITRANGE)
430  | QPMS_COORDS_POL;
431  return res;
432  }
433 }
434 
436 static inline PGenReturnDataBulk PGen_fetch_cart2(struct PGen *g,
437  size_t nmemb, cart2_t *target) {
438  if (g->c->fetch_cart2)
439  return g->c->fetch_cart2(g, nmemb, target);
440  else {
441  anycoord_point_t *tmp;
442  QPMS_CRASHING_MALLOC(tmp, sizeof(anycoord_point_t) * nmemb);
443  PGenReturnDataBulk res = PGen_fetch_any(g, nmemb, tmp);
445  tmp, res.flags, res.generated);
446  free(tmp);
447  res.flags = (res.flags & ~QPMS_COORDS_BITRANGE)
449  return res;
450  }
451 }
452 
454 static inline PGenReturnDataBulk PGen_fetch_z(struct PGen *g,
455  size_t nmemb, double *target) {
456  if (g->c->fetch_z)
457  return g->c->fetch_z(g, nmemb, target);
458  else {
459  anycoord_point_t *tmp;
460  QPMS_CRASHING_MALLOC(tmp, sizeof(anycoord_point_t) * nmemb);
461  PGenReturnDataBulk res = PGen_fetch_any(g, nmemb, tmp);
463  tmp, res.flags, res.generated);
464  free(tmp);
465  res.flags = (res.flags & ~QPMS_COORDS_BITRANGE)
467  return res;
468  }
469 }
470 
472 static inline void PGen_destroy(PGen *g) {
473  g->c->destructor(g);
474  assert(g->stateData == NULL); // this should be done by the destructor
475 }
476 
478 static inline PGenZReturnData PGen_next_z(PGen *g) {
479  if (g->c->next_z)
480  return g->c->next_z(g);
481  else { // Super-slow generic fallback.
482  PGenReturnData res = PGen_next(g);
483  if (res.flags & PGEN_NOTDONE) {
484  PGenZReturnData r;
485  r.point_z = anycoord2cart1(res.point, res.flags);
486  r.flags = (res.flags & ~QPMS_COORDS_BITRANGE)
488  return r;
489  } else
490  return PGenZDoneVal;
491  }
492 }
493 
496  if (g->c->next_sph)
497  return g->c->next_sph(g);
498  else { // Super-slow generic fallback.
499  PGenReturnData res = PGen_next(g);
500  if (res.flags & PGEN_NOTDONE) {
502  r.point_sph = anycoord2sph(res.point, res.flags);
503  r.flags = (res.flags & ~QPMS_COORDS_BITRANGE)
504  | QPMS_COORDS_SPH;
505  return r;
506  } else
507  return PGenSphDoneVal;
508  }
509 }
510 
513  if (g->c->next_pol)
514  return g->c->next_pol(g);
515  else { // Super-slow generic fallback.
516  PGenReturnData res = PGen_next(g);
517  if (res.flags & PGEN_NOTDONE) {
519  r.point_pol = anycoord2pol(res.point, res.flags);
520  r.flags = (res.flags & ~QPMS_COORDS_BITRANGE)
521  | QPMS_COORDS_POL;
522  return r;
523  } else
524  return PGenPolDoneVal;
525  }
526 }
527 
530  if (g->c->next_cart3)
531  return g->c->next_cart3(g);
532  else { // Super-slow generic fallback.
533  PGenReturnData res = PGen_next(g);
534  if (res.flags & PGEN_NOTDONE) {
536  r.point_cart3 = anycoord2cart3(res.point, res.flags);
537  r.flags = (res.flags & ~QPMS_COORDS_BITRANGE)
539  return r;
540  } else
541  return PGenCart3DoneVal;
542  }
543 }
544 
547  if (g->c->next_cart2)
548  return g->c->next_cart2(g);
549  else { // Super-slow generic fallback.
550  PGenReturnData res = PGen_next(g);
551  if (res.flags & PGEN_NOTDONE) {
553  r.point_cart2 = anycoord2cart2(res.point, res.flags);
554  r.flags = (res.flags & ~QPMS_COORDS_BITRANGE)
556  return r;
557  } else
558  return PGenCart2DoneVal;
559  }
560 }
561 
562 #if 0
564 static inline PGenReturnDataBulk PGen_fetch(PGen *g, size_t n, anycoord_point_t *arr){
565  if (g->c->fetch)
566  return g->c->fetch(g, n, arr);
567  else abort();
568 }
569 #endif
570 
571 static inline bool PGenSph_notDone(PGenSphReturnData data) {
572  return data.flags & PGEN_NOTDONE ? true : false;
573 }
574 static inline bool PGenCart3_notDone(PGenCart3ReturnData data) {
575  return data.flags & PGEN_NOTDONE ? true : false;
576 }
577 
580  PGenCart3ReturnData c3data;
581  c3data.flags = sphdata.flags;
582  c3data.point_cart3 = sph2cart(sphdata.point_sph);
583  return c3data;
584 }
585 
589  s.flags = c.flags;
590  s.point_sph = cart2sph(c.point_cart3);
591  return s;
592 }
593 
594 /*
595  * Some basic lattice generators implementing the abstract interface above (implemented in latticegens.c).
596  */
597 
599 extern const PGenClassInfo PGen_FromPoint2DArray; // TODO Do I even need this to be declared here?
601 PGen PGen_FromPoints2DArray_new(const point2d *points, size_t len);
602 
604 extern const PGenClassInfo PGen_1D;
605 typedef enum PGen_1D_incrementDirection{
606  //PGEN_1D_POSITIVE_INC, // not implemented
607  //PGEN_1D_NEGATIVE_INC, // not implemented
608  PGEN_1D_INC_FROM_ORIGIN,
609  PGEN_1D_INC_TOWARDS_ORIGIN
610 } PGen_1D_incrementDirection;
612 PGen PGen_1D_new_minMaxR(double period,
613  double offset,
614  double minR,
615  bool inc_minR,
616  double maxR,
617  bool inc_maxR,
618  PGen_1D_incrementDirection incdir
619 );
620 
621 extern const PGenClassInfo PGen_xyWeb;
622 PGen PGen_xyWeb_new(cart2_t b1, cart2_t b2, double rtol, cart2_t offset,
623  double minR, bool inc_minR, double maxR, bool inc_maxR);
624 
626 size_t PGen_xyWeb_sizecap(cart2_t b1, cart2_t b2, double rtol, cart2_t offset,
627  double minR, bool inc_minR, double maxR, bool inc_maxR);
628 
629 
630 extern const PGenClassInfo PGen_LatticeRadialHeap2D;
631 extern const PGenClassInfo PGen_LatticeRadialHeap3D;
632 PGen PGen_LatticeRadialHeap2D_new(cart2_t b1, cart2_t b2, cart2_t offset,
633  double minR, bool inc_minR, double maxR, bool inc_maxR);
634 PGen PGen_LatticeRadialHeap3D_new(const cart3_t *b1, const cart3_t *b2, const cart3_t *b3,
635  const cart3_t *offset, double minR, bool inc_minR, double maxR, bool inc_maxR);
636 
637 
639 extern const PGenClassInfo PGen_shifted;
640 PGen PGen_shifted_new(PGen orig, cart3_t shift);
641 
642 /*
643  * THE NICE PART (adaptation of lattices2d.py)
644  * ===========================================
645  *
646  * all the functions are prefixed with l2d_
647  * convention for argument order: inputs, *outputs, add. params
648  */
649 
650 #define BASIS_RTOL 1e-13
651 
652 // Bravais lattice types
653 typedef enum {
654  OBLIQUE = 1,
655  RECTANGULAR = 2,
656  SQUARE = 4,
657  RHOMBIC = 5,
658  EQUILATERAL_TRIANGULAR = 3,
659  RIGHT_ISOSCELES=SQUARE,
660  PARALLELOGRAMMIC=OBLIQUE,
661  CENTERED_RHOMBIC=RECTANGULAR,
662  RIGHT_TRIANGULAR=RECTANGULAR,
663  CENTERED_RECTANGULAR=RHOMBIC,
664  ISOSCELE_TRIANGULAR=RHOMBIC,
665  RIGHT_ISOSCELE_TRIANGULAR=SQUARE,
666  HEXAGONAL=EQUILATERAL_TRIANGULAR
667 } LatticeType2;
668 
669 #if 0
670 // Wallpaper groups
671 typedef enum {
672  TODO
673 } SpaceGroup2;
674 #endif
675 
676 // Just for detecting the right angles (needed for generators).
677 typedef enum {
678  NOT_ORTHOGONAL = 0,
679  ORTHOGONAL_01 = 1,
680  ORTHOGONAL_12 = 2,
681  ORTHOGONAL_02 = 4
682 } LatticeFlags;
683 
684 
685 
686 /*
687  * Lagrange-Gauss reduction of a 2D basis.
688  * The output shall satisfy |out1| <= |out2| <= |out2 - out1|
689  */
690 void l2d_reduceBasis(cart2_t in1, cart2_t in2, cart2_t *out1, cart2_t *out2);
691 
692 // This one uses LLL reduction.
693 void l3d_reduceBasis(const cart3_t in[3], cart3_t out[3]);
694 
695 /*
696  * This gives the "ordered shortest triple" of base vectors (each pair from the triple
697  * is a base) and there may not be obtuse angle between o1, o2 and between o2, o3
698  */
699 void l2d_shortestBase3(cart2_t i1, cart2_t i2, cart2_t *o1, cart2_t *o2, cart2_t *o3);
700 
701 /*
702  * TODO doc
703  * return value is 4 or 6.
704  */
705 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);
706 // variant
707 int l2d_shortestBase46_arr(cart2_t i1, cart2_t i2, cart2_t *oarr, double rtol);
708 
709 // Determines whether angle between inputs is obtuse
710 bool l2d_is_obtuse_r(cart2_t i1, cart2_t i2, double rtol);
711 bool l2d_is_obtuse(cart2_t i1, cart2_t i2);
712 
713 /*
714  * Given two basis vectors, returns 2D Bravais lattice type.
715  */
716 LatticeType2 l2d_classifyLattice(cart2_t b1, cart2_t b2, double rtol);
717 
718 // Detects right angles.
719 LatticeFlags l2d_detectRightAngles(cart2_t b1, cart2_t b2, double rtol);
720 LatticeFlags l3d_detectRightAngles(const cart3_t basis[3], double rtol);
721 
722 // Other functions in lattices2d.py: TODO?
723 // range2D()
724 // generateLattice()
725 // generateLatticeDisk()
726 // cutWS()
727 // filledWS()
728 // change_basis()
729 
730 /*
731  * Given basis vectors, returns the corners of the Wigner-Seits unit cell (W1, W2, -W1, W2)
732  * for rectangular and square lattice or (w1, w2, w3, -w1, -w2, -w3) otherwise.
733  */
734 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);
735 // variant
736 int l2d_cellCornersWS_arr(cart2_t i1, cart2_t i2, cart2_t *oarr, double rtol);
737 
738 // Reciprocal bases; returns 0 on success, possibly a non-zero if b1 and b2 are parallel
739 int l2d_reciprocalBasis1(cart2_t b1, cart2_t b2, cart2_t *rb1, cart2_t *rb2);
740 int l2d_reciprocalBasis2pi(cart2_t b1, cart2_t b2, cart2_t *rb1, cart2_t *rb2);
741 // 3D reciprocal bases; returns (direct) unit cell volume with possible sign. Assumes direct lattice basis already reduced.
742 double l3d_reciprocalBasis1(const cart3_t direct_basis[3], cart3_t reciprocal_basis[3]);
743 double l3d_reciprocalBasis2pi(const cart3_t direct_basis[3], cart3_t reciprocal_basis[3]);
744 
745 double l2d_unitcell_area(cart2_t b1, cart2_t b2);
746 
747 // returns the radius of inscribed circle of a hexagon (or rectangle/square if applicable) created by the shortest base triple
748 double l2d_hexWebInCircleRadius(cart2_t b1, cart2_t b2);
749 
750 /*
751  * THE MORE OR LESS OK PART
752  * ========================
753  */
754 
755 /*
756  * General set of points ordered by the r-coordinate.
757  * Typically, this will include all lattice inside a certain circle.
758  * This structure is internally used by the "lattice generators" below.
759  * It does not have its memory management of its own, as it is handled
760  * by the "generators". For everything except the generators,
761  * this structure shall be read-only.
762  */
763 typedef struct {
764  size_t nrs; // number of different radii
765  double *rs; // the radii; of length nrs (largest contained radius == rs[nrs-1])
766  point2d *base;
767  ptrdiff_t *r_offsets; // of length nrs+1 (using relative offsets due to possible realloc's)
768  // the jth point of i-th radius is base[r_offsets[i]+j] or using the inline below..
769  /* // redundand (therefore removed) members
770  * point2d *points; // redundant as it is the same as points_at_r[0]
771  * size_t npoints; // redundant as it is the same as points_at_r[nrs]-points_at_r[0]
772  */
774 
775 
776 // sorts arbitrary points and creates points2d_rordered_t
777 points2d_rordered_t *points2d_rordered_frompoints(const point2d *orig_base,
778  size_t nmemb, double rtol, double atol);
779 
780 // returns a copy but shifted by a constant (actually in a stupid way, but whatever)
781 points2d_rordered_t *points2d_rordered_shift(const points2d_rordered_t *orig,
782  point2d shift, double rtol, double atol);
783 
784 // returns a copy but scaled by a factor
785 points2d_rordered_t *points2d_rordered_scale(const points2d_rordered_t *orig,
786  double factor);
787 
788 
789 /* The destructor: use only for results of
790  * - points2D_rordered_frompoints,
791  * - points2d_rordered_shift,
792  * - points2d_rordered_scale.
793  */
794 void points2d_rordered_free(points2d_rordered_t *);
795 
796 static inline point2d points2d_rordered_get_point(const points2d_rordered_t *ps, int r_order, int i) {
797  assert(i >= 0);
798  assert(r_order < ps->nrs);
799  assert(i < (ps->r_offsets[r_order+1] - ps->r_offsets[r_order]));
800  return ps->base[ps->r_offsets[r_order] + i];
801 }
802 
803 static inline double points2d_rordered_get_r(const points2d_rordered_t *ps, int r_order) {
804  assert(r_order < ps->nrs);
805  return ps->rs[r_order];
806 }
807 
808 
809 ptrdiff_t points2d_rordered_locate_r(const points2d_rordered_t *, double r);
810 
811 // returns a "view" (does not copy any of the arrays)
812 // -- DO NOT FREE orig BEFORE THE END OF SCOPE OF THE RESULT
813 points2d_rordered_t points2d_rordered_annulus(const points2d_rordered_t *orig, double minr, bool minr_inc,
814  double maxr, bool maxr_inc);
815 
816 
817 #if 0
819 double qpms_emptylattice2_mode_nth(
820  cart2_t b1_rec,
821  cart2_t b2_rec,
822  double rtol,
823  cart2_t k,
824  double wave_speed,
825  size_t N
826  );
827 
829 void qpms_emptylattice2_modes_maxindex(
830  double target_freqs[],
831  cart2_t b1_rec,
832  cart2_t b2_rec,
833  double rtol,
834  cart2_t k,
835  double wave_speed,
836  size_t maxindex
837  );
838 #endif
839 
841 
848  double **target_freqs,
849  cart2_t b1_rec,
850  cart2_t b2_rec,
851  double rtol,
852  cart2_t k,
853  double wave_speed,
854  double maxfreq
855  );
856 
859  double target[2],
860  cart2_t b1_rec,
861  cart2_t b2_rec,
862  double rtol,
863  cart2_t k,
864  double wave_speed,
865  double omega
866  );
867 
868 
869 
870 /*
871  * THE UGLY PART
872  * =============
873  */
874 
875 /*
876  * EQUILATERAL TRIANGULAR LATTICE
877  */
878 
879 typedef enum {
880  TRIANGULAR_VERTICAL, // there is a lattice base vector parallel to the y-axis
881  TRIANGULAR_HORIZONTAL // there is a lattice base vector parallel to the x-axis
882 } TriangularLatticeOrientation;
883 
884 static inline TriangularLatticeOrientation reverseTriangularLatticeOrientation(TriangularLatticeOrientation o){
885  switch(o) {
886  case TRIANGULAR_VERTICAL:
887  return TRIANGULAR_HORIZONTAL;
888  break;
889  case TRIANGULAR_HORIZONTAL:
890  return TRIANGULAR_VERTICAL;
891  break;
892  default:
893  abort();
894  }
895  abort();
896 }
897 
898 // implementation data structures; not needed in the header file
900 
901 typedef struct {
902  // public:
904  TriangularLatticeOrientation orientation;
905  double a; // lattice vector length
906 
907  // not sure if needed:
908  bool includes_origin;
909 
910  // Denotes an offset of the "origin" point; meaning step hexshift * a / sqrt(2) upwards
911  // or leftwards for the horizontal or vertical orientations, respectively.
912  int hexshift;
913 
914  // private:
916 
918 
919 triangular_lattice_gen_t *triangular_lattice_gen_init(double a, TriangularLatticeOrientation ori, bool include_origin,
920  int halfoffset);
921 const points2d_rordered_t * triangular_lattice_gen_getpoints(const triangular_lattice_gen_t *g);
922 int triangular_lattice_gen_extend_to_r(triangular_lattice_gen_t *g, double r);
923 int triangular_lattice_gen_extend_to_steps(triangular_lattice_gen_t *g, int maxsteps);
924 void triangular_lattice_gen_free(triangular_lattice_gen_t *g);
925 
926 
927 /*
928  * HONEYCOMB LATTICE
929  */
930 
931 typedef struct {
932  // public:
934  TriangularLatticeOrientation orientation;
935  double a;
936  double h;
937 
938  // private:
941 
942 honeycomb_lattice_gen_t *honeycomb_lattice_gen_init_h(double h, TriangularLatticeOrientation ori);
943 honeycomb_lattice_gen_t *honeycomb_lattice_gen_init_a(double a, TriangularLatticeOrientation ori);
944 int honeycomb_lattice_gen_extend_to_steps(honeycomb_lattice_gen_t *g, int maxsteps);
945 int honeycomb_lattice_gen_extend_to_r(honeycomb_lattice_gen_t *g, double r);
946 void honeycomb_lattice_gen_free(honeycomb_lattice_gen_t *g);
947 
948 #endif // LATTICES_H
const PGenClassInfo PGen_FromPoint2DArray
2D point generator that simply iterates over an existing array of Point2d.
Definition: latticegens.c:234
static PGenReturnData PGen_next_nf(struct PGen *g)
Generate a point with any of the next-methods.
Definition: lattices.h:268
static PGenCart3ReturnData PGen_next_cart3(PGen *g)
Generate a point in a 3D real space (cartesian coordinates).
Definition: lattices.h:529
PGen PGen_1D_new_minMaxR(double period, double offset, double minR, bool inc_minR, double maxR, bool inc_maxR, PGen_1D_incrementDirection incdir)
PGen_1D point generator constructor.
Definition: latticegens.c:296
static PGenReturnDataBulk PGen_fetch_z(struct PGen *g, size_t nmemb, double *target)
Generate multiple points in 1D cartesian coordinates.
Definition: lattices.h:454
PGenPointFlags
Definition: lattices.h:127
@ PGEN_METHOD_UNAVAILABLE
Set if no suitable method exists (no point generated).
Definition: lattices.h:149
@ PGEN_DONE
Convenience identifier, not an actual flag.
Definition: lattices.h:150
@ PGEN_SINGLE_R
Definition: lattices.h:146
@ PGEN_NOTDONE
Definition: lattices.h:133
@ PGEN_AT_Z
Set if the point(s) lie(s) at the z-axis (theta is either 0 or M_PI).
Definition: lattices.h:147
@ PGEN_AT_XY
Set if the point(s) lie(s) in the xy-plane (theta is M_PI2).
Definition: lattices.h:148
@ PGEN_OLD_R
Definition: lattices.h:139
struct PGenReturnDataBulk PGenReturnDataBulk
Metadata generated by the fetch*() methods from PGenClassInfo.
static PGenReturnDataBulk PGen_fetch_sph(struct PGen *g, size_t nmemb, sph_t *target)
Generate multiple points in spherical coordinates.
Definition: lattices.h:382
static PGenReturnDataBulk PGen_fetch_cart2(struct PGen *g, size_t nmemb, cart2_t *target)
Generate multiple points in 2D cartesian coordinates.
Definition: lattices.h:436
struct PGenClassInfo PGenClassInfo
PGen class metadata.
struct PGenReturnData PGenReturnData
Generic PGen return type that might contain point represented in any of the supported coordinate syst...
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)
Gives the frequencies of two empty lattice modes nearest to omega at a given wave vector k.
Definition: lattices2d.c:929
static PGenCart3ReturnData PGenReturnDataConv_sph_cart3(PGenSphReturnData sphdata)
Standard PGen spherical coordinates -> 3d cartesian convertor.
Definition: lattices.h:579
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)
Gives the frequencies of empty lattice modes at a given wave vector k up to maxfreq and one more.
Definition: lattices2d.c:886
struct PGen PGen
Generic lattice point generator type.
struct PGenSphReturnData PGenSphReturnData
PGen single-point return data type (3D, spherical coordinates).
int qpms_reduce_lattice_basis(double *b, const size_t bsize, const size_t ndim, double delta)
Lattice basis reduction.
Definition: lll.c:41
static PGenSphReturnData PGen_next_sph(PGen *g)
Generate a point in a 3D real space (spherical coordinates).
Definition: lattices.h:495
static void PGen_destroy(PGen *g)
Deallocate and invalidate a PGen point generator.
Definition: lattices.h:472
struct PGenPolReturnData PGenPolReturnData
PGen single-point return data type (2D, polar coordinates).
struct PGenCart2ReturnData PGenCart2ReturnData
PGen single-point return data type (2D, cartesian coordinates).
static PGenPolReturnData PGen_next_pol(PGen *g)
Generate a point in a 2D real space (polar coordinates).
Definition: lattices.h:512
static PGenReturnDataBulk PGen_fetch_cart3(struct PGen *g, size_t nmemb, cart3_t *target)
Generate multiple points in 3D cartesian coordinates.
Definition: lattices.h:400
struct PGenZReturnData PGenZReturnData
PGen single-point return data type (1D).
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)
Returns a number larger or equal than the number of all the points generated by a PGen_xyWeb.
Definition: latticegens.c:591
static PGenReturnDataBulk PGen_fetch_pol(struct PGen *g, size_t nmemb, pol_t *target)
Generate multiple points in polar coordinates.
Definition: lattices.h:418
const PGenClassInfo PGen_1D
1D equidistant point generator.
Definition: latticegens.c:407
const PGenClassInfo PGen_shifted
A metagenerator generating points from another generator shifted by a constant.
Definition: latticegens.c:992
static PGenZReturnData PGen_next_z(PGen *g)
Generate a point in a 1D real space.
Definition: lattices.h:478
static PGenCart2ReturnData PGen_next_cart2(PGen *g)
Generate a point in a 2D real space (cartesian coordinates).
Definition: lattices.h:546
PGen PGen_FromPoints2DArray_new(const point2d *points, size_t len)
PGen_FromPoint2DArray constructor.
static PGenSphReturnData PGenReturnDataConv_cart3_sph(PGenCart3ReturnData c)
Standard PGen 3d cartesian -> spherical coordinates convertor.
Definition: lattices.h:587
static PGenReturnData PGen_next(struct PGen *g)
Generate a point with any of the next-methods or fetch-methods.
Definition: lattices.h:366
struct PGenCart3ReturnData PGenCart3ReturnData
PGen single-point return data type (3D, cartesian coordinates).
static PGenReturnDataBulk PGen_fetch_any(struct PGen *g, size_t nmemb, anycoord_point_t *target)
Generate multiple points with PGen in any coordinate system.
Definition: lattices.h:301
#define QPMS_CRASHING_MALLOC(pointer, size)
Wrapper macro of standard malloc(), crashing on failure.
Definition: qpms_error.h:109
#define QPMS_ENSURE(x, msg,...)
Raises an error if x is not true, with custom error message.
Definition: qpms_error.h:226
@ QPMS_COORDS_CART2
2D cartesian.
Definition: qpms_types.h:244
@ QPMS_COORDS_CART3
Definition: qpms_types.h:245
@ QPMS_COORDS_BITRANGE
Convenience bitmask (not a valid flag!).
Definition: qpms_types.h:247
@ QPMS_COORDS_CART1
1D cartesian (= double).
Definition: qpms_types.h:241
@ QPMS_COORDS_SPH
3D spherical.
Definition: qpms_types.h:243
@ QPMS_COORDS_POL
2D polar.
Definition: qpms_types.h:242
union anycoord_point_t anycoord_point_t
Union type capable to contain various 1D, 2D and 3D coordinates.
PGen single-point return data type (2D, cartesian coordinates).
Definition: lattices.h:192
cart2_t point_cart2
Generated point in cartesian coordinates.
Definition: lattices.h:194
PGenPointFlags flags
Metadata.
Definition: lattices.h:193
PGen single-point return data type (3D, cartesian coordinates).
Definition: lattices.h:198
PGenPointFlags flags
Metadata.
Definition: lattices.h:199
cart3_t point_cart3
Generated point in cartesian coordinates.
Definition: lattices.h:200
PGen class metadata.
Definition: lattices.h:231
PGenReturnDataBulk(* fetch_cart3)(struct PGen *, size_t, cart3_t *)
Generate up to n 3D points in cartesian coordinates.
Definition: lattices.h:259
PGenReturnDataBulk(* fetch)(struct PGen *, size_t, anycoord_point_t *)
Generate up to n points in the native coordinates.
Definition: lattices.h:249
void(* destructor)(struct PGen *)
Destructor.
Definition: lattices.h:264
PGenPointFlags native_point_flags
Info about the generator native coordinate system.
Definition: lattices.h:235
PGenReturnDataBulk(* fetch_z)(struct PGen *, size_t, double *)
Generate up to n 1D points.
Definition: lattices.h:251
PGenReturnDataBulk(* fetch_cart2)(struct PGen *, size_t, cart2_t *)
Generate up to n 2D points in cartesian coordinates.
Definition: lattices.h:257
PGenPolReturnData(* next_pol)(struct PGen *)
Generate a single 2D point in polar coordinates.
Definition: lattices.h:241
PGenCart2ReturnData(* next_cart2)(struct PGen *)
Generate a single 2D point in cartesian coordinates.
Definition: lattices.h:245
PGenSphReturnData(* next_sph)(struct PGen *)
Generate a single 3D point in spherical coordinates.
Definition: lattices.h:243
PGenReturnDataBulk(* fetch_pol)(struct PGen *, size_t, pol_t *)
Generate up to n 2D points in polar coordinates.
Definition: lattices.h:253
PGenZReturnData(* next_z)(struct PGen *)
Generate a single 1D point.
Definition: lattices.h:239
PGenCart3ReturnData(* next_cart3)(struct PGen *)
Generate a single 3D point in cartesian coordinates.
Definition: lattices.h:247
PGenReturnDataBulk(* fetch_sph)(struct PGen *, size_t, sph_t *)
Generate up to n 3D points in spherical coordinates.
Definition: lattices.h:255
PGenReturnData(* next)(struct PGen *)
Generate a single point in the native coordinates.
Definition: lattices.h:237
PGen single-point return data type (2D, polar coordinates).
Definition: lattices.h:180
pol_t point_pol
Generated point in polar coordinates.
Definition: lattices.h:182
PGenPointFlags flags
Metadata.
Definition: lattices.h:181
Metadata generated by the fetch*() methods from PGenClassInfo.
Definition: lattices.h:161
size_t generated
Number of points really generated.
Definition: lattices.h:164
PGenPointFlags flags
Flags describing the returned data.
Definition: lattices.h:163
Generic PGen return type that might contain point represented in any of the supported coordinate syst...
Definition: lattices.h:168
anycoord_point_t point
Generated point in a coordinate system defined by flags.
Definition: lattices.h:170
PGenPointFlags flags
Metadata, must contain valid coordinate system defining flags.
Definition: lattices.h:169
PGen single-point return data type (3D, spherical coordinates).
Definition: lattices.h:186
sph_t point_sph
Generated point in spherical coordinates.
Definition: lattices.h:188
PGenPointFlags flags
Metadata.
Definition: lattices.h:187
PGen single-point return data type (1D).
Definition: lattices.h:174
double point_z
Generated point on a real axis.
Definition: lattices.h:176
PGenPointFlags flags
Medatata.
Definition: lattices.h:175
Generic lattice point generator type.
Definition: lattices.h:120
void * stateData
Pointer to internal state data; shall be NULL if invalid (destroyed);.
Definition: lattices.h:124
const struct PGenClassInfo * c
Pointer to the "class" metadata defining the behaviour of the generator.
Definition: lattices.h:122
2D cartesian coordinates. See also vectors.h.
Definition: qpms_types.h:202
3D cartesian coordinates. See also vectors.h.
Definition: qpms_types.h:186
Definition: lattices.h:931
Definition: lattices.h:763
2D polar coordinates. See also vectors.h.
Definition: qpms_types.h:223
Spherical coordinates. See also vectors.h.
Definition: qpms_types.h:207
Definition: lattices2d.c:258
Definition: lattices.h:901
Union type capable to contain various 1D, 2D and 3D coordinates.
Definition: qpms_types.h:229
double z
1D cartesian coordinate.
Definition: qpms_types.h:230
Coordinate transforms and vector arithmetics.
static pol_t anycoord2pol(anycoord_point_t p, qpms_coord_system_t t)
Conversion from anycoord_point_t to explicitly polar coordinates.
Definition: vectors.h:523
static sph_t cart2sph(const cart3_t cart)
3D cartesian to spherical coordinates conversion. See Coordinate systems and default conversions.
Definition: vectors.h:134
static cart3_t anycoord2cart3(anycoord_point_t p, qpms_coord_system_t t)
Conversion from anycoord_point_t to explicitly 3D cartesian coordinates.
Definition: vectors.h:484
static void anycoord_arr2something(void *dest, qpms_coord_system_t tdest, const anycoord_point_t *src, qpms_coord_system_t tsrc, size_t nmemb)
Coordinate conversion of point arrays (anycoord_point_t to something).
Definition: vectors.h:746
static sph_t anycoord2sph(anycoord_point_t p, qpms_coord_system_t t)
Conversion from anycoord_point_t to explicitly spherical coordinates.
Definition: vectors.h:459
static cart3_t sph2cart(const sph_t sph)
Spherical to 3D cartesian coordinates conversion. See Coordinate systems and default conversions.
Definition: vectors.h:143
static double anycoord2cart1(anycoord_point_t p, qpms_coord_system_t t)
Conversion from anycoord_point_t to explicitly 1D cartesian coordinates.
Definition: vectors.h:573
static cart2_t anycoord2cart2(anycoord_point_t p, qpms_coord_system_t t)
Conversion from anycoord_point_t to explicitly 2D cartesian coordinates.
Definition: vectors.h:548