Experimental support for periodic lattices in scatsys "constructor".

Former-commit-id: 69727f4d415866b948af83c55ed8adb46b651f16

qpms/lattices.h

@@ -726,6 +726,9 @@ 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);
 

qpms/lattices2d.c

@@ -849,6 +849,21 @@ int l2d_reciprocalBasis2pi(cart2_t b1, cart2_t b2, cart2_t *rb1, cart2_t *rb2) {
   return retval;
 };
 
+// 3D reciprocal bases; returns (direct) unit cell volume. Assumes direct lattice basis already reduced.
+double l3d_reciprocalBasis1(const cart3_t db[3], cart3_t rb[3]) {
+  double vol = cart3_tripleprod(db[0], db[1], db[2]);
+  for(int i = 0; i < 3; ++i)
+    rb[i] = cart3_divscale(cart3_vectorprod(db[(i+1) % 3], db[(i+2) % 3]), vol);
+  return vol;
+}
+
+double l3d_reciprocalBasis2pi(const cart3_t db[3], cart3_t rb[3]) {
+  double vol = l3d_reciprocalBasis1(db, rb);
+  for(int i = 1; i < 3; ++i)
+    rb[i] = cart3_scale(2 * M_PI, rb[i]);
+  return vol;
+}
+
 // 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 i1, cart2_t i2) {
   cart2_t b1, b2, b3;

qpms/scatsystem.c

@@ -142,15 +142,68 @@ static void add_orbit_type(qpms_scatsys_t *ss, const qpms_ss_orbit_type_t *ot_cu
   ss->orbit_type_count++;
 }
 
+// Standardises the lattice base vectors, and fills the other contents of ss->per[0].
+// LPTODO split the functionality in smaller functions, these might be useful elsewhere.
+static void process_lattice_bases(qpms_scatsys_t *ss, const qpms_tolerance_spec_t *tol) {
+  switch(ss->lattice_dimension) {
+    case 1: {
+              double normsq;
+              normsq = cart3_normsq(ss->per->lattice_basis[0]);
+              ss->per->unitcell_volume = sqrt(normsq);
+              ss->per->reciprocal_basis[0] = cart3_divscale(ss->per->lattice_basis[0], normsq);
+      } break;
+    case 2: {
+              // (I) Create orthonormal basis of the plane in which the basis vectors lie
+              cart3_t onbasis[2];
+              double norm0 = cart3norm(ss->per->lattice_basis[0]);
+              // 0: Just normalised 0. basis vector
+              onbasis[0] = cart3_divscale(ss->per->lattice_basis[0], norm0); 
+              // 1: Gram-Schmidt complement of the other basis vector
+              const double b0norm_dot_b1 = cart3_dot(onbasis[0], ss->per->lattice_basis[1]);
+              onbasis[1] = cart3_substract(ss->per->lattice_basis[1],
+                  cart3_scale(b0norm_dot_b1, onbasis[0]));
+              onbasis[1] = cart3_divscale(onbasis[1], cart3norm(onbasis[1]));
+              // (II) Express the lattice basis in terms of the new 2d plane vector basis onbasis
+              cart2_t b2d[2] = {{.x = norm0, .y = 0}, 
+                {.x = b0norm_dot_b1, .y = cart3_dot(onbasis[1], ss->per->lattice_basis[1])}};
+              // (III) Reduce lattice basis and get reciprocal bases in the 2D plane
+              l2d_reduceBasis(b2d[0], b2d[1], &b2d[0], &b2d[1]);
+              ss->per->unitcell_volume = l2d_unitcell_area(b2d[0], b2d[1]);
+              cart2_t r2d[2];
+              QPMS_ENSURE_SUCCESS(l2d_reciprocalBasis1(b2d[0], b2d[1], &r2d[0], &r2d[1]));
+              // (IV) Rotate everything back to the original 3D space.
+              for(int i = 0; i < 2; ++i) {
+                ss->per->lattice_basis[i] = cart3_add(
+                    cart3_scale(b2d[i].x, onbasis[0]),
+                    cart3_scale(b2d[i].y, onbasis[1]));
+                ss->per->reciprocal_basis[i] = cart3_add(
+                    cart3_scale(r2d[i].x, onbasis[0]),
+                    cart3_scale(r2d[i].y, onbasis[1]));
+              }
+      } break;
+    case 3: {
+              l3d_reduceBasis(ss->per->lattice_basis, ss->per->lattice_basis); 
+              ss->per->unitcell_volume = fabs(l3d_reciprocalBasis1(ss->per->lattice_basis,
+                  ss->per->reciprocal_basis));
+              // TODO check unitcell_volume for sanity w.r.t tolerance
+      } break;
+    default:
+      QPMS_WTF;
+  }
+}
+
 
 // Almost 200 lines. This whole thing deserves a rewrite!
 qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig,
     const qpms_finite_group_t *sym, complex double omega, 
     const qpms_tolerance_spec_t *tol) {
-  if(orig->lattice_dimension)
-    QPMS_NOT_IMPLEMENTED("Periodic systems not yet done");
-  // TODO check data sanity
+  if (sym == NULL) {
+    QPMS_WARN("No symmetry group pointer provided, assuming trivial symmetry"
+        " (this is currently the only option for periodic systems).");
+    sym = &QPMS_FINITE_GROUP_TRIVIAL;
+  }
 
+  // TODO check data sanity
 
   qpms_l_t lMax = 0; // the overall lMax of all base specs.
   qpms_normalisation_t normalisation = QPMS_NORMALISATION_UNDEF;
@@ -180,9 +233,15 @@ qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig,
 
   // Allocate T-matrix, particle and particle orbit info arrays
   qpms_scatsys_t *ss;
-  QPMS_CRASHING_MALLOC(ss, sizeof(qpms_scatsys_t));     
+  QPMS_CRASHING_MALLOC(ss, sizeof(*ss) + !!orig->lattice_dimension * sizeof(ss->per[0]));
   ss->lattice_dimension = orig->lattice_dimension;
   // TODO basic periodic lattices related stuff here.
+  if (ss->lattice_dimension) {
+    ss->per[0] = orig->per[0];
+    process_lattice_bases(ss, tol);
+  }
+  for(int i = 0; i < ss->lattice_dimension; ++i)  // extend lenscale by basis vectors
+    lenscale = MAX(lenscale, cart3norm(ss->per->lattice_basis[i]));
   ss->lenscale = lenscale;
   ss->sym = sym;
   ss->medium = orig->medium;

qpms/scatsystem.h

@@ -132,18 +132,13 @@ typedef struct qpms_ss_derived_tmatrix_t {
 } qpms_ss_derived_tmatrix_t;
 
 typedef struct qpms_scatsys_periodic_info_t {	
-	/// Coordinate system for \a lattice_basis.
-	/** This is mandatory for \a lattice_dimension != 0 */
-	qpms_coord_system_t lattice_basis_csystem;
 	/// (Direct) lattice basis of the system (only \a lattice_dimension elements are used)
 	/** This is mandatory for \a lattice_dimension != 0 */
-	anycoord_point_t lattice_basis[3];
+	cart3_t lattice_basis[3];
 
-	/// Coordinate system for \a reciprocal_basis.
-	qpms_coord_system_t reciprocal_basis_csystem;
 	/// Reciprocal lattice basis. 
 	/**(TODO specify whether it includes 2π or not) */
-	anycoord_point_t reciprocal_basis[3];
+	cart3_t reciprocal_basis[3];
 
 	/// Unitcell volume (irrelevant for non-periodic systems).
 	/** The dimensionality of the volume corresponds to lattice_dimension, so
@@ -151,7 +146,7 @@ typedef struct qpms_scatsys_periodic_info_t {
 	 * lattice_dimension == 2, a 2D area.
 	 */
 	double unitcell_volume;
-} qpms_scatsys_pediodic_info_t;
+} qpms_scatsys_periodic_info_t;
 
 
 struct qpms_trans_calculator;
@@ -160,14 +155,8 @@ struct qpms_epsmu_generator_t;
 /// Common "class" for system of scatterers, both periodic and non-periodic.
 /**
  * Infinite periodic structures (those with \a lattice_dimension > 0) 
- * have the following additional members filled:
- *  - lattice_basis_csystem
- *  - lattice_basis
- *  - reciprocal_basis_csystem
- *  - reciprocal_basis
- *  - unitcell_volume
+ * have the \a per element allocated and filled.
  * These are ignored for finite systems (lattice_dimension == 0).
- *
  */
 typedef struct qpms_scatsys_t {
 	/// Number of dimensions in which the system is periodic from the range 0–3.
@@ -225,7 +214,7 @@ typedef struct qpms_scatsys_t {
 	double lenscale; // radius of the array, used as a relative tolerance measure
 	struct qpms_trans_calculator *c;
 
-	/// Periodic lattice metadata. Only allocated/used when lattice_dimension != 0.
+	/// Periodic lattice metadata. Only allocated/used when lattice_dimension != 0 (exactly one member).
 	qpms_scatsys_periodic_info_t per[];
 } qpms_scatsys_t;
 
@@ -256,6 +245,7 @@ typedef struct qpms_scatsys_at_omega_t {
  * so keep them alive until scatsys is destroyed.
  *  
  *  The following fields must be filled in the "proto- scattering system" \a orig:
+ *  * orig->lattice_dimension
  *  * orig->medium – The pointers are copied to the new qpms_scatsys_t instance; 
  *    the target qpms_abstract_tmatrix_t objects must be kept alive before all the resulting 
  *    qpms_scatsys_t instances are properly destroyed.
@@ -270,6 +260,12 @@ typedef struct qpms_scatsys_at_omega_t {
  *  * orig->p
  *  * orig->p_count
  *
+ *  For periodic systems, the corresponding number of orig->per->lattice_basis[] elements
+ *  must be filled as well.
+ *
+ *  For periodic systems, only trivial group is currently supported. Non-trivial 
+ *  groups will cause undefined behaviour.
+ *
  *  The resulting qpms_scatsys_t is obtained by actually evaluating the T-matrices
  *  at the given frequency \a omega and where applicable, these are compared
  *  by their values with given tolerances. The T-matrix generators are expected

qpms/vectors.h

@@ -104,6 +104,22 @@ static inline double cart3_dot(const cart3_t a, const cart3_t b) {
 	return a.x * b.x + a.y * b.y + a.z * b.z;
 }
 
+
+/// 3D vector product a.k.a. cross product.
+static inline cart3_t cart3_vectorprod(const cart3_t a, const cart3_t b) {
+	cart3_t c = {
+		.x = a.y * b.z - a.z * b.y,
+		.y = a.z * b.x - a.x * b.z,
+		.z = a.x * b.y - a.y * b.x,
+	};
+	return c;
+}
+
+/// Scalar triple product \f$ a \cdot ( b \times c ) \f$.
+static inline double cart3_tripleprod(const cart3_t a, const cart3_t b, const cart3_t c) {
+	return cart3_dot(a, cart3_vectorprod(b, c));
+}
+
 /// 3D vector euclidian norm squared.
 static inline double cart3_normsq(const cart3_t a) {
 	return cart3_dot(a, a);
@@ -171,6 +187,12 @@ static inline cart3_t cart3_scale(const double c, const cart3_t v) {
 	return res;
 }
 
+/// 3D vector division by scalar (N.B. argument order).
+static inline cart3_t cart3_divscale( const cart3_t v, const double c) {
+	cart3_t res = {v.x / c, v.y / c, v.z / c};
+	return res;
+}
+
 /// Euclidian distance between two 3D points.
 static inline double cart3_dist(const cart3_t a, const cart3_t b) {
 	return cart3norm(cart3_substract(a,b));
@@ -480,6 +502,12 @@ static inline cart3_t anycoord2cart3(anycoord_point_t p, qpms_coord_system_t t)
 	QPMS_WTF;
 }
 
+/// Cartesian norm of anycoord_point_t.
+// The implementation is simple and stupid, do not use for heavy computations.
+static inline double anycoord_norm(anycoord_point_t p, qpms_coord_system_t t) {
+	return cart3norm(anycoord2cart3(p, t));
+}
+
 #if 0
 // Convenience identifiers for return values.
 static const cart3_t CART3_INVALID = {NAN, NAN, NAN};
@@ -846,7 +874,16 @@ static inline void anycoord_arr2something(void *dest, qpms_coord_system_t tdest,
   }
 }
 
+/// Converts cart3_t to array of doubles.
+static inline void cart3_to_double_array(double a[], cart3_t b) {
+	a[0] = b.x; a[1] = b.y; a[2] = b.z;
+}
 
+/// Converts array of doubles to cart3_t.
+static inline cart3_t cart3_from_double_array(const double a[]) {
+	cart3_t b = {.x = a[0], .y = a[1], .z = a[1]};
+	return b;
+}
 
 
 typedef double matrix3d[3][3];