From 0f03509dde6fd059cec8ca110bee0fdfcabe2103 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 13 Sep 2019 11:52:56 +0300 Subject: [PATCH] Mode frequency calculations for empty lattices. Former-commit-id: ca14f8d765399c5bf92bb64cc7f1dc6722d53ebb --- qpms/latticegens.c | 15 +++++++++ qpms/lattices.h | 59 ++++++++++++++++++++++++++++++++-- qpms/lattices2d.c | 79 ++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 151 insertions(+), 2 deletions(-) diff --git a/qpms/latticegens.c b/qpms/latticegens.c index 6eed602..44dfd8f 100644 --- a/qpms/latticegens.c +++ b/qpms/latticegens.c @@ -585,3 +585,18 @@ const PGenClassInfo PGen_xyWeb = { }; +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) +{ + l2d_reduceBasis(b1, b2, &b1, &b2); + LatticeFlags lf = l2d_detectRightAngles(b1, b2, rtol); + double layer_min_height = l2d_hexWebInCircleRadius(b1, b2); + long layer = ceil(minR / layer_min_height); + if(!inc_minR && (layer * layer_min_height) <= minR) + ++layer; + long last_layer = floor(maxR / layer_min_height); + if(!inc_maxR && (last_layer * layer_min_height) >= maxR) + --(last_layer); + // TODO less crude estimate (this one should be safe, however) + return ((lf & ORTHOGONAL_01) ? 4 : 6) * (last_layer - layer + 1); +} diff --git a/qpms/lattices.h b/qpms/lattices.h index bbf515a..01ce9d7 100644 --- a/qpms/lattices.h +++ b/qpms/lattices.h @@ -542,7 +542,7 @@ static inline PGenCart3ReturnData PGen_next_cart3(PGen *g) { } } -/// Ǧenerate a point in a 2D real space (cartesian coordinates). +/// Generate a point in a 2D real space (cartesian coordinates). static inline PGenCart2ReturnData PGen_next_cart2(PGen *g) { if (g->c->next_cart2) return g->c->next_cart2(g); @@ -618,11 +618,14 @@ PGen PGen_1D_new_minMaxR(double period, ///< Distance between points. PGen_1D_incrementDirection incdir ///< Order of generated points. ); - 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); +/// Returns a number larger or equal than the number of all the points generated by a PGen_xyWeb. +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); + /* * THE NICE PART (adaptation of lattices2d.py) @@ -796,6 +799,58 @@ points2d_rordered_t points2d_rordered_annulus(const points2d_rordered_t *orig, d double maxr, bool maxr_inc); +#if 0 +/// Gives the frequency of \a n-th empty lattice mode at a given wave vector \a k. +double qpms_emptylattice2_mode_nth( + cart2_t b1_rec, ///< First reciprocal lattice base vector + cart2_t b2_rec, ///< Second reciprocal lattice base vector + double rtol, ///< Relative tolerance to detect right angles + cart2_t k, ///< The wave vector + double wave_speed, ///< Wave speed in a given medium (i.e. vacuum speed / refractive index). + size_t N ///< Index of the mode (note that degenerate modes are counted multiple times). + ); + +/// Gives the first `maxindex` frequencies of empty lattice modes at a given wave vector \a k. +void qpms_emptylattice2_modes_maxindex( + double target_freqs[], ///< Target array of size maxindex. + cart2_t b1_rec, ///< First reciprocal lattice base vector + cart2_t b2_rec, ///< Second reciprocal lattice base vector + double rtol, ///< Relative tolerance to detect right angles + cart2_t k, ///< The wave vector + double wave_speed, ///< Wave speed in a given medium (i.e. vacuum speed / refractive index). + size_t maxindex ///< Number of the frequencies generated. + ); +#endif + +/// Gives the frequencies of empty lattice modes at a given wave vector \a k up to \a maxfreq and one more. +/** + * The frequencies are saved to a newly allocated array *target_freqs (to be deallocated + * using free() by the caller). + * + * \returns Number of found mode frequencies lower or equal than \a maxfreq plus one. + */ +size_t qpms_emptylattice2_modes_maxfreq( + double **target_freqs, + cart2_t b1_rec, ///< First reciprocal lattice base vector + cart2_t b2_rec, ///< Second reciprocal lattice base vector + double rtol, ///< Relative tolerance to detect right angles + cart2_t k, ///< The wave vector + double wave_speed, ///< Wave speed in a given medium (i.e. vacuum speed / refractive index). + double maxfreq ///< The maximum frequency. + ); + +/// Gives the frequencies of two empty lattice modes nearest to \a omega at a given wave vector \a k. +void qpms_emptylattice2_modes_nearest( + double target[2], ///< Target array with lower ([0]) and upper ([1]) frequency. + cart2_t b1_rec, ///< First reciprocal lattice base vector + cart2_t b2_rec, ///< Second reciprocal lattice base vector + double rtol, ///< Relative tolerance to detect right angles + cart2_t k, ///< The wave vector + double wave_speed, ///< Wave speed in a given medium (i.e. vacuum speed / refractive index). + double omega ///< The frequency around which the frequencies are searched. + ); + + /* * THE UGLY PART diff --git a/qpms/lattices2d.c b/qpms/lattices2d.c index 3a9ef69..4b94929 100644 --- a/qpms/lattices2d.c +++ b/qpms/lattices2d.c @@ -866,4 +866,83 @@ double l2d_unitcell_area(cart2_t b1, cart2_t b2) { } +#if 0 +double qpms_emptylattice2_mode_nth( + cart2_t b1_rec, + cart2_t b2_rec, + double rtol, + cart2_t k, + double c, + 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 c, + size_t maxindex + ); +#endif + +static int dblcmp(const void *p1, const void *p2) { + const double *x1 = (double *) p1, *x2 = (double *) p2; + double dif = *x1 - *x2; + if(dif > 0) return 1; + else if (dif < 0) return -1; + else return 0; +} + +size_t qpms_emptylattice2_modes_maxfreq(double **target_freqs, + cart2_t b1, cart2_t b2, double rtol, cart2_t k, + double c, double maxfreq) +{ + QPMS_UNTESTED; + l2d_reduceBasis(b1, b2, &b1, &b2); + double maxk_safe = cart2norm(k) + maxfreq / c + cart2norm(b1) + cart2norm(b2); // This is an overkill + size_t capacity = PGen_xyWeb_sizecap(b1, b2, rtol, k, 0, true, maxk_safe, true); + cart2_t *Kpoints; + QPMS_CRASHING_MALLOC(Kpoints, sizeof(cart2_t) * capacity); + PGen Kgen = PGen_xyWeb_new(b1, b2, rtol, k, 0, true, maxk_safe, true); + + size_t generated; + PGenReturnDataBulk rd; + while((rd = PGen_fetch_cart2(&Kgen, capacity - generated, Kpoints + generated)).flags & PGEN_NOTDONE) + generated += rd.generated; + + double *thefreqs; + QPMS_CRASHING_MALLOC(thefreqs, generated * sizeof(double)); + for(size_t i = 0; i < generated; ++i) thefreqs[i] = cart2norm(Kpoints[i]) * c; + + free(Kpoints); + + qsort(thefreqs, generated, sizeof(double), dblcmp); + size_t count; + bool hitmax = false; + for(count = 0; count < generated; ++count) + if(thefreqs[count] > maxfreq) { + if(hitmax) + break; + else + hitmax = true; + } + + *target_freqs = thefreqs; + return count; +} + +void qpms_emptylattice2_modes_nearest(double target[2], + cart2_t b1, cart2_t b2, double rtol, + cart2_t k, double c, double omega) +{ + QPMS_UNTESTED; + double *freqlist; + size_t n = qpms_emptylattice2_modes_maxfreq(&freqlist, + b1, b2, rtol, k, c, omega); + target[0] = freqlist[n-2]; + target[1] = freqlist[n-1]; + free(freqlist); +}