diff --git a/qpms/lattices.h b/qpms/lattices.h index 01ce9d7..aa5fc5e 100644 --- a/qpms/lattices.h +++ b/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); diff --git a/qpms/lattices2d.c b/qpms/lattices2d.c index 6249f00..edb2f47 100644 --- a/qpms/lattices2d.c +++ b/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; diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index e4965d0..8b4e7fa 100644 --- a/qpms/scatsystem.c +++ b/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; diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 3d34aa1..051d9c4 100644 --- a/qpms/scatsystem.h +++ b/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 diff --git a/qpms/vectors.h b/qpms/vectors.h index cf4b01c..679aaf3 100644 --- a/qpms/vectors.h +++ b/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];