2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2009,2010,2011,2012,2013 by the GROMACS development team.
5 * Copyright (c) 2014,2015,2016,2017,2018 by the GROMACS development team.
6 * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
38 * \page page_module_selection_insolidangle Selection method: insolidangle
40 * This method selects a subset of particles that are located in a solid
41 * angle defined by a center and a set of points.
42 * The solid angle is constructed as a union of small cones whose axis
43 * goes through the center and a point.
44 * So there's such a cone for each position, and a
45 * point is in the solid angle if it lies within any of these cones.
46 * The width of the cones can be adjusted.
48 * The method is implemented by partitioning the surface of the unit sphere
49 * into bins using the polar coordinates \f$(\theta, \phi)\f$.
50 * The partitioning is always uniform in the zenith angle \f$\theta\f$,
51 * while the partitioning in the azimuthal angle \f$\phi\f$ varies.
52 * For each reference point, the unit vector from the center to the point
53 * is constructed, and it is stored in all the bins that overlap with the
54 * cone defined by the point.
55 * Bins that are completely covered by a single cone are marked as such.
56 * Checking whether a point is in the solid angle is then straightforward
57 * with this data structure: one finds the bin that corresponds to the point,
58 * and checks whether the bin is completely covered. If it is not, one
59 * additionally needs to check whether it is within the specified cutoff of
60 * any of the stored points.
62 * The above construction gives quite a lot of flexibility for constructing
63 * the bins without modifying the rest of the code.
64 * The current (quite inefficient) implementation is discussed below, but
65 * it should be optimized to get the most out of the code.
67 * The current way of constructing the bins constructs the boundaries
68 * statically: the bin size in the zenith direction is set to approximately
69 * half the angle cutoff, and the bins in the azimuthal direction have
70 * sizes such that the shortest edge of the bin is approximately equal to
71 * half the angle cutoff (for the regions close to the poles, a single bin
73 * Each reference point is then added to the bins as follows:
74 * -# Find the zenith angle range that is spanned by the cone centered at the
75 * point (this is simple addition/subtraction).
76 * -# Calculate the maximal span of the cone in the azimuthal direction using
78 * \f[\sin \Delta \phi_{max} = \frac{\sin \alpha}{\sin \theta}\f]
79 * (a sine formula in spherical coordinates),
80 * where \f$\alpha\f$ is the width of the cone and \f$\theta\f$ is the
81 * zenith angle of the cone center.
82 * Similarly, the zenith angle at which this extent is achieved is
84 * \f[\cos \theta_{max} = \frac{\cos \theta}{\cos \alpha}\f]
85 * (Pythagoras's theorem in spherical coordinates).
86 * -# For each zenith angle bin that is at least partially covered by the
87 * cone, calculate the span of the cone at the edges using
88 * \f[\sin^2 \frac{\Delta \phi}{2} = \frac{\sin^2 \frac{\alpha}{2} - \sin^2 \frac{\theta -
89 * \theta'}{2}}{\sin \theta \sin \theta'}\f] (distance in spherical geometry), where \f$\theta'\f$
90 * is the zenith angle of the bin edge. Treat zenith angle bins that are completely covered by the
91 * cone (in the case that the cone is centered close to the pole) as a special case.
92 * -# Using the values calculated above, loop through the azimuthal bins that
93 * are partially or completely covered by the cone and update them.
95 * The total solid angle (for covered fraction calculations) is estimated by
96 * taking the total area of completely covered bins plus
97 * half the area of partially covered bins.
98 * The second one is an approximation, but should give reasonable estimates
99 * for the averages as well as in cases where the bin size is small.
103 * Implements the \ref sm_insolidangle "insolidangle" selection method.
106 * The implementation could be optimized quite a bit.
109 * Move the covered fraction stuff somewhere else and make it more generic
110 * (along the lines it is handled in selection.h and trajana.h in the old C
113 * \author Teemu Murtola <teemu.murtola@gmail.com>
114 * \ingroup module_selection
122 #include "gromacs/math/functions.h"
123 #include "gromacs/math/units.h"
124 #include "gromacs/math/utilities.h"
125 #include "gromacs/math/vec.h"
126 #include "gromacs/pbcutil/pbc.h"
127 #include "gromacs/selection/indexutil.h"
128 #include "gromacs/selection/selection.h"
129 #include "gromacs/utility/arraysize.h"
130 #include "gromacs/utility/exceptions.h"
131 #include "gromacs/utility/smalloc.h"
133 #include "position.h"
135 #include "selmethod.h"
136 #include "selmethod_impl.h"
143 * Internal data structure for the \p insolidangle selection method.
145 * \see \c t_partition
147 * \ingroup module_selection
151 /** Left edge of the partition. */
153 /** Bin index corresponding to this partition. */
159 * Internal data structure for the \p insolidangle selection method.
161 * Describes the surface partitioning within one slice along the zenith angle.
162 * The slice from azimuthal angle \p p[i].left to \p p[i+1].left belongs to
165 * \ingroup module_selection
167 typedef struct partition
169 /** Number of partition items (\p p contains \p n+1 items). */
171 /** Array of partition edges and corresponding bins. */
177 * Internal data structure for the \p insolidangle selection method.
179 * Contains the reference points that partially cover a certain region on the
180 * surface of the unit sphere.
181 * If \p n is -1, the whole region described by the bin is covered.
183 * \ingroup module_selection
185 typedef struct spheresurfacebin
187 /** Number of points in the array \p x, -1 if whole bin covered. */
189 /** Number of elements allocated for \p x. */
191 /** Array of points that partially cover the bin. */
193 } t_spheresurfacebin;
197 * Data structure for the \p insolidangle selection method.
199 * All angle values are in the units of radians.
201 * \ingroup module_selection
203 typedef struct methoddata_insolidangle
205 /** Center of the solid angle. */
206 gmx_ana_pos_t center;
207 /** Positions that span the solid angle. */
211 /** Estimate of the covered fraction. */
214 /** Cutoff for the cosine (equals cos(angcut)). */
216 /** Bin size to be used as the target bin size when constructing the bins. */
219 /** Number of bins in the \p tbin array. */
221 /** Size of one bin in the zenith angle direction. */
223 /** Array of zenith angle slices. */
225 /** Number of elements allocated for the \p bin array. */
227 /** Number of elements used in the \p bin array. */
229 /** Array of individual bins. */
230 t_spheresurfacebin* bin;
231 } t_methoddata_insolidangle;
234 * Allocates data for the \p insolidangle selection method.
236 * \param[in] npar Not used (should be 3).
237 * \param[in,out] param Method parameters (should point to
238 * \ref smparams_insolidangle).
239 * \returns Pointer to the allocated data (\ref t_methoddata_insolidangle).
241 * Allocates memory for a \ref t_methoddata_insolidangle structure and
242 * initializes the parameter as follows:
243 * - \p center defines the value for t_methoddata_insolidangle::center.
244 * - \p span defines the value for t_methoddata_insolidangle::span.
245 * - \p cutoff defines the value for t_methoddata_insolidangle::angcut.
247 static void* init_data_insolidangle(int npar, gmx_ana_selparam_t* param);
249 * Initializes the \p insolidangle selection method.
251 * \param top Not used.
252 * \param npar Not used.
253 * \param param Not used.
254 * \param data Pointer to \ref t_methoddata_insolidangle to initialize.
255 * \returns 0 on success, -1 on failure.
257 * Converts t_methoddata_insolidangle::angcut to radians and allocates
258 * and allocates memory for the bins used during the evaluation.
260 static void init_insolidangle(const gmx_mtop_t* top, int npar, gmx_ana_selparam_t* param, void* data);
261 /** Frees the data allocated for the \p insolidangle selection method. */
262 static void free_data_insolidangle(void* data);
264 * Initializes the evaluation of the \p insolidangle selection method for a frame.
266 * \param[in] context Evaluation context.
267 * \param data Should point to a \ref t_methoddata_insolidangle.
269 * Creates a lookup structure that enables fast queries of whether a point
270 * is within the solid angle or not.
272 static void init_frame_insolidangle(const gmx::SelMethodEvalContext& context, void* data);
273 /** Internal helper function for evaluate_insolidangle(). */
274 static bool accept_insolidangle(rvec x, const t_pbc* pbc, void* data);
275 /** Evaluates the \p insolidangle selection method. */
276 static void evaluate_insolidangle(const gmx::SelMethodEvalContext& context,
278 gmx_ana_selvalue_t* out,
281 /** Calculates the distance between unit vectors. */
282 static real sph_distc(rvec x1, rvec x2);
283 /** Does a binary search on a \p t_partition to find a bin for a value. */
284 static int find_partition_bin(t_partition* p, real value);
285 /** Finds a bin that corresponds to a location on the unit sphere surface. */
286 static int find_surface_bin(t_methoddata_insolidangle* surf, rvec x);
287 /** Clears/initializes the bins on the unit sphere surface. */
288 static void clear_surface_points(t_methoddata_insolidangle* surf);
289 /** Frees memory allocated for storing the reference points in the surface bins. */
290 static void free_surface_points(t_methoddata_insolidangle* surf);
291 /** Adds a reference point to a given bin. */
292 static void add_surface_point(t_methoddata_insolidangle* surf, int tbin, int pbin, rvec x);
293 /** Marks a bin as completely covered. */
294 static void mark_surface_covered(t_methoddata_insolidangle* surf, int tbin, int pbin);
295 /** Helper function for store_surface_point() to update a single zenith angle bin. */
296 static void update_surface_bin(t_methoddata_insolidangle* surf,
303 /** Adds a single reference point and updates the surface bins. */
304 static void store_surface_point(t_methoddata_insolidangle* surf, rvec x);
306 * Optimizes the surface bins for faster searching.
308 * \param[in,out] surf Surface data structure.
310 * Currently, this function does nothing.
312 static void optimize_surface_points(t_methoddata_insolidangle* surf);
313 /** Estimates the area covered by the reference cones. */
314 static real estimate_covered_fraction(t_methoddata_insolidangle* surf);
315 /** Checks whether a point lies within a solid angle. */
316 static bool is_surface_covered(t_methoddata_insolidangle* surf, rvec x);
318 /** Parameters for the \p insolidangle selection method. */
319 static gmx_ana_selparam_t smparams_insolidangle[] = {
320 { "center", { POS_VALUE, 1, { nullptr } }, nullptr, SPAR_DYNAMIC },
321 { "span", { POS_VALUE, -1, { nullptr } }, nullptr, SPAR_DYNAMIC | SPAR_VARNUM },
322 { "cutoff", { REAL_VALUE, 1, { nullptr } }, nullptr, SPAR_OPTIONAL },
325 /** Help text for the \p insolidangle selection method. */
326 static const char* const help_insolidangle[] = {
329 " insolidangle center POS span POS_EXPR [cutoff REAL]",
331 "This keyword selects atoms that are within [TT]REAL[tt] degrees",
332 "(default=5) of any position in [TT]POS_EXPR[tt] as seen from [TT]POS[tt]",
333 "a position expression that evaluates to a single position), i.e., atoms",
334 "in the solid angle spanned by the positions in [TT]POS_EXPR[tt] and",
335 "centered at [TT]POS[tt].[PAR]",
337 "Technically, the solid angle is constructed as a union of small cones",
338 "whose tip is at [TT]POS[tt] and the axis goes through a point in",
339 "[TT]POS_EXPR[tt]. There is such a cone for each position in",
340 "[TT]POS_EXPR[tt], and point is in the solid angle if it lies within any",
341 "of these cones. The cutoff determines the width of the cones.",
344 /** Selection method data for the \p insolidangle method. */
345 gmx_ana_selmethod_t sm_insolidangle = {
349 asize(smparams_insolidangle),
350 smparams_insolidangle,
351 &init_data_insolidangle,
355 &free_data_insolidangle,
356 &init_frame_insolidangle,
358 &evaluate_insolidangle,
359 { "insolidangle center POS span POS_EXPR [cutoff REAL]",
360 "Selecting atoms in a solid angle",
361 asize(help_insolidangle),
365 static void* init_data_insolidangle(int /* npar */, gmx_ana_selparam_t* param)
367 t_methoddata_insolidangle* data = new t_methoddata_insolidangle();
371 data->distccut = 0.0;
372 data->targetbinsize = 0.0;
375 data->tbinsize = 0.0;
376 data->tbin = nullptr;
381 param[0].val.u.p = &data->center;
382 param[1].val.u.p = &data->span;
383 param[2].val.u.r = &data->angcut;
387 static void init_insolidangle(const gmx_mtop_t* /* top */,
389 gmx_ana_selparam_t* /* param */,
392 t_methoddata_insolidangle* surf = static_cast<t_methoddata_insolidangle*>(data);
395 if (surf->angcut <= 0)
397 GMX_THROW(gmx::InvalidInputError("Angle cutoff should be > 0"));
400 surf->angcut *= gmx::c_deg2Rad;
402 surf->distccut = -std::cos(surf->angcut);
403 surf->targetbinsize = surf->angcut / 2;
404 surf->ntbins = static_cast<int>(M_PI / surf->targetbinsize);
405 surf->tbinsize = (180.0 / surf->ntbins) * gmx::c_deg2Rad;
407 snew(surf->tbin, static_cast<int>(M_PI / surf->tbinsize) + 1);
409 for (i = 0; i < surf->ntbins; ++i)
411 c = static_cast<int>(std::max(std::sin(surf->tbinsize * i), std::sin(surf->tbinsize * (i + 1)))
412 * M_2PI / surf->targetbinsize)
414 snew(surf->tbin[i].p, c + 1);
418 snew(surf->bin, surf->maxbins);
422 * \param data Data to free (should point to a \ref t_methoddata_insolidangle).
424 * Frees the memory allocated for \c t_methoddata_insolidangle::center and
425 * \c t_methoddata_insolidangle::span, as well as the memory for the internal
428 static void free_data_insolidangle(void* data)
430 t_methoddata_insolidangle* d = static_cast<t_methoddata_insolidangle*>(data);
435 for (i = 0; i < d->ntbins; ++i)
441 free_surface_points(d);
446 static void init_frame_insolidangle(const gmx::SelMethodEvalContext& context, void* data)
448 t_methoddata_insolidangle* d = static_cast<t_methoddata_insolidangle*>(data);
452 free_surface_points(d);
453 clear_surface_points(d);
454 for (i = 0; i < d->span.count(); ++i)
458 pbc_dx(context.pbc, d->span.x[i], d->center.x[0], dx);
462 rvec_sub(d->span.x[i], d->center.x[0], dx);
465 store_surface_point(d, dx);
467 optimize_surface_points(d);
472 * \param[in] x Test point.
473 * \param[in] pbc PBC data (if NULL, no PBC are used).
474 * \param[in] data Pointer to a \c t_methoddata_insolidangle data structure.
475 * \returns true if \p x is within the solid angle, false otherwise.
477 static bool accept_insolidangle(rvec x, const t_pbc* pbc, void* data)
479 t_methoddata_insolidangle* d = static_cast<t_methoddata_insolidangle*>(data);
484 pbc_dx(pbc, x, d->center.x[0], dx);
488 rvec_sub(x, d->center.x[0], dx);
491 return is_surface_covered(d, dx);
495 * See sel_updatefunc() for description of the parameters.
496 * \p data should point to a \c t_methoddata_insolidangle.
498 * Calculates which atoms in \p g are within the solid angle spanned by
499 * \c t_methoddata_insolidangle::span and centered at
500 * \c t_methoddata_insolidangle::center, and stores the result in \p out->u.g.
502 static void evaluate_insolidangle(const gmx::SelMethodEvalContext& context,
504 gmx_ana_selvalue_t* out,
508 for (int b = 0; b < pos->count(); ++b)
510 if (accept_insolidangle(pos->x[b], context.pbc, data))
512 gmx_ana_pos_add_to_group(out->u.g, pos, b);
518 * \param[in] sel Selection element to query.
519 * \returns true if the covered fraction can be estimated for \p sel with
520 * _gmx_selelem_estimate_coverfrac(), false otherwise.
522 bool _gmx_selelem_can_estimate_cover(const gmx::SelectionTreeElement& sel)
524 if (sel.type == SEL_BOOLEAN && sel.u.boolt == BOOL_OR)
529 bool bDynFound = false;
530 gmx::SelectionTreeElementPointer child = sel.child;
533 if (child->type == SEL_EXPRESSION)
535 if (child->u.expr.method->name == sm_insolidangle.name)
537 if (bFound || bDynFound)
543 else if (child->u.expr.method && (child->u.expr.method->flags & SMETH_DYNAMIC))
552 else if (!_gmx_selelem_can_estimate_cover(*child))
562 * \param[in] sel Selection for which the fraction should be calculated.
563 * \returns Fraction of angles covered by the selection (between zero and one).
565 * The return value is undefined if _gmx_selelem_can_estimate_cover() returns
567 * Should be called after gmx_ana_evaluate_selections() has been called for the
570 real _gmx_selelem_estimate_coverfrac(const gmx::SelectionTreeElement& sel)
574 if (sel.type == SEL_EXPRESSION && sel.u.expr.method->name == sm_insolidangle.name)
576 t_methoddata_insolidangle* d = static_cast<t_methoddata_insolidangle*>(sel.u.expr.mdata);
579 d->cfrac = estimate_covered_fraction(d);
583 if (sel.type == SEL_BOOLEAN && sel.u.boolt == BOOL_NOT)
585 cfrac = _gmx_selelem_estimate_coverfrac(*sel.child);
593 /* Here, we assume that the selection is simple enough */
594 gmx::SelectionTreeElementPointer child = sel.child;
597 cfrac = _gmx_selelem_estimate_coverfrac(*child);
608 * \param[in] x1 Unit vector 1.
609 * \param[in] x2 Unit vector 2.
610 * \returns Minus the dot product of \p x1 and \p x2.
612 * This function is used internally to calculate the distance between the
613 * unit vectors \p x1 and \p x2 to find out whether \p x2 is within the
614 * cone centered at \p x1. Currently, the cosine of the angle is used
615 * for efficiency, and the minus is there to make it behave like a normal
616 * distance (larger values mean longer distances).
618 static real sph_distc(rvec x1, rvec x2)
620 return -iprod(x1, x2);
624 * \param[in] p Partition to search.
625 * \param[in] value Value to search for.
626 * \returns The partition index in \p p that contains \p value.
628 * If \p value is outside the range of \p p, the first/last index is returned.
629 * Otherwise, the return value \c i satisfies \c p->p[i].left<=value and
630 * \c p->p[i+1].left>value
632 static int find_partition_bin(t_partition* p, real value)
634 int pmin, pmax, pbin;
636 /* Binary search the partition */
639 while (pmax > pmin + 1)
641 pbin = pmin + (pmax - pmin) / 2;
642 if (p->p[pbin].left <= value)
656 * \param[in] surf Surface data structure to search.
657 * \param[in] x Unit vector to find.
658 * \returns The bin index that contains \p x.
660 * The return value is an index to the \p surf->bin array.
662 static int find_surface_bin(t_methoddata_insolidangle* surf, rvec x)
668 phi = atan2(x[YY], x[XX]);
669 tbin = static_cast<int>(std::floor(theta / surf->tbinsize));
670 if (tbin >= surf->ntbins)
672 tbin = surf->ntbins - 1;
674 pbin = find_partition_bin(&surf->tbin[tbin], phi);
675 return surf->tbin[tbin].p[pbin].bin;
679 * \param[in,out] surf Surface data structure.
681 * Clears the reference points from the bins and (re)initializes the edges
682 * of the azimuthal bins.
684 static void clear_surface_points(t_methoddata_insolidangle* surf)
689 for (i = 0; i < surf->ntbins; ++i)
691 c = static_cast<int>(std::min(std::sin(surf->tbinsize * i), std::sin(surf->tbinsize * (i + 1)))
692 * M_2PI / surf->targetbinsize)
699 for (j = 0; j < c; ++j)
701 surf->tbin[i].p[j].left = -M_PI + j * M_2PI / c - 0.0001;
702 surf->tbin[i].p[j].bin = surf->nbins;
703 surf->bin[surf->nbins].n = 0;
706 surf->tbin[i].p[c].left = M_PI + 0.0001;
707 surf->tbin[i].p[c].bin = -1;
712 * \param[in,out] surf Surface data structure.
714 static void free_surface_points(t_methoddata_insolidangle* surf)
718 for (i = 0; i < surf->nbins; ++i)
722 sfree(surf->bin[i].x);
724 surf->bin[i].n_alloc = 0;
725 surf->bin[i].x = nullptr;
730 * \param[in,out] surf Surface data structure.
731 * \param[in] tbin Bin number in the zenith angle direction.
732 * \param[in] pbin Bin number in the azimuthal angle direction.
733 * \param[in] x Point to store.
735 static void add_surface_point(t_methoddata_insolidangle* surf, int tbin, int pbin, rvec x)
739 bin = surf->tbin[tbin].p[pbin].bin;
740 /* Return if bin is already completely covered */
741 if (surf->bin[bin].n == -1)
745 /* Allocate more space if necessary */
746 if (surf->bin[bin].n == surf->bin[bin].n_alloc)
748 surf->bin[bin].n_alloc += 10;
749 srenew(surf->bin[bin].x, surf->bin[bin].n_alloc);
751 /* Add the point to the bin */
752 copy_rvec(x, surf->bin[bin].x[surf->bin[bin].n]);
757 * \param[in,out] surf Surface data structure.
758 * \param[in] tbin Bin number in the zenith angle direction.
759 * \param[in] pbin Bin number in the azimuthal angle direction.
761 static void mark_surface_covered(t_methoddata_insolidangle* surf, int tbin, int pbin)
765 bin = surf->tbin[tbin].p[pbin].bin;
766 surf->bin[bin].n = -1;
770 * \param[in,out] surf Surface data structure.
771 * \param[in] tbin Bin number in the zenith angle direction.
772 * \param[in] phi Azimuthal angle of \p x.
773 * \param[in] pdelta1 Width of the cone at the lower edge of \p tbin.
774 * \param[in] pdelta2 Width of the cone at the uppper edge of \p tbin.
775 * \param[in] pdeltamax Max. width of the cone inside \p tbin.
776 * \param[in] x Point to store (should have unit length).
778 static void update_surface_bin(t_methoddata_insolidangle* surf,
786 real pdelta, phi1, phi2;
787 int pbin1, pbin2, pbiniter, pbin;
789 /* Find the edges of the bins affected */
790 pdelta = max(max(pdelta1, pdelta2), pdeltamax);
794 pbin = find_partition_bin(&surf->tbin[tbin], phi1);
799 pbin = find_partition_bin(&surf->tbin[tbin], phi1 + M_2PI);
800 pbin1 = pbin - surf->tbin[tbin].n;
805 pbin2 = find_partition_bin(&surf->tbin[tbin], phi2);
809 pbin2 = find_partition_bin(&surf->tbin[tbin], phi2 - M_2PI);
810 pbin2 += surf->tbin[tbin].n;
813 if (pbin2 - pbin1 > surf->tbin[tbin].n)
815 pbin2 = pbin1 + surf->tbin[tbin].n;
817 /* Find the edges of completely covered region */
818 pdelta = min(pdelta1, pdelta2);
825 /* Loop over all affected bins */
826 for (pbiniter = pbin1; pbiniter != pbin2; ++pbiniter, ++pbin)
828 /* Wrap bin around if end reached */
829 if (pbin == surf->tbin[tbin].n)
835 /* Check if bin is completely covered and update */
836 if (surf->tbin[tbin].p[pbin].left >= phi1 && surf->tbin[tbin].p[pbin + 1].left <= phi2)
838 mark_surface_covered(surf, tbin, pbin);
842 add_surface_point(surf, tbin, pbin, x);
848 * \param[in,out] surf Surface data structure.
849 * \param[in] x Point to store (should have unit length).
851 * Finds all the bins covered by the cone centered at \p x and calls
852 * update_surface_bin() to update them.
854 static void store_surface_point(t_methoddata_insolidangle* surf, rvec x)
857 real pdeltamax, tmax;
858 real theta1, theta2, pdelta1, pdelta2;
862 phi = atan2(x[YY], x[XX]);
863 /* Find the maximum extent in the phi direction */
864 if (theta <= surf->angcut)
869 else if (theta >= M_PI - surf->angcut)
876 pdeltamax = std::asin(sin(surf->angcut) / sin(theta));
877 tmax = std::acos(cos(theta) / cos(surf->angcut));
879 /* Find the first affected bin */
880 tbin = max(static_cast<int>(std::floor((theta - surf->angcut) / surf->tbinsize)), 0);
881 theta1 = tbin * surf->tbinsize;
882 if (theta1 < theta - surf->angcut)
890 /* Loop through all affected bins */
891 while (tbin < std::ceil((theta + surf->angcut) / surf->tbinsize) && tbin < surf->ntbins)
893 /* Calculate the next boundaries */
894 theta2 = (tbin + 1) * surf->tbinsize;
895 if (theta2 > theta + surf->angcut)
897 /* The circle is completely outside the cone */
900 else if (theta2 <= -(theta - surf->angcut) || theta2 >= M_2PI - (theta + surf->angcut)
901 || tbin == surf->ntbins - 1)
903 /* The circle is completely inside the cone, or we are in the
904 * 360 degree bin covering the pole. */
909 /* TODO: This formula is numerically unstable if theta is very
910 * close to the pole. In practice, it probably does not matter
911 * much, but it would be nicer to adjust the theta bin boundaries
912 * such that the case above catches this instead of falling through
915 * asin(std::sqrt((gmx::square(std::sin(surf->angcut / 2))
916 - gmx::square(std::sin((theta2 - theta) / 2)))
917 / (sin(theta) * sin(theta2))));
920 if (tmax >= theta1 && tmax <= theta2)
922 update_surface_bin(surf, tbin, phi, pdelta1, pdelta2, pdeltamax, x);
926 update_surface_bin(surf, tbin, phi, pdelta1, pdelta2, 0, x);
935 static void optimize_surface_points(t_methoddata_insolidangle* /* surf */)
937 /* TODO: Implement */
941 * \param[in] surf Surface data structure.
942 * \returns An estimate for the area covered by the reference points.
944 static real estimate_covered_fraction(t_methoddata_insolidangle* surf)
947 real cfrac, tfrac, pfrac;
950 for (t = 0; t < surf->ntbins; ++t)
952 tfrac = std::cos(t * surf->tbinsize) - std::cos((t + 1) * surf->tbinsize);
953 for (p = 0; p < surf->tbin[t].n; ++p)
955 pfrac = surf->tbin[t].p[p + 1].left - surf->tbin[t].p[p].left;
956 n = surf->bin[surf->tbin[t].p[p].bin].n;
957 if (n == -1) /* Bin completely covered */
959 cfrac += tfrac * pfrac;
961 else if (n > 0) /* Bin partially covered */
963 cfrac += tfrac * pfrac / 2; /* A rough estimate */
967 return cfrac / (4 * M_PI);
971 * \param[in] surf Surface data structure to search.
972 * \param[in] x Unit vector to check.
973 * \returns true if \p x is within the solid angle, false otherwise.
975 static bool is_surface_covered(t_methoddata_insolidangle* surf, rvec x)
979 bin = find_surface_bin(surf, x);
980 /* Check for completely covered bin */
981 if (surf->bin[bin].n == -1)
985 /* Check each point that partially covers the bin */
986 for (i = 0; i < surf->bin[bin].n; ++i)
988 if (sph_distc(x, surf->bin[bin].x[i]) < surf->distccut)