2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2009,2010,2011,2012,2013,2014,2015, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
36 * \page page_module_selection_insolidangle Selection method: insolidangle
38 * This method selects a subset of particles that are located in a solid
39 * angle defined by a center and a set of points.
40 * The solid angle is constructed as a union of small cones whose axis
41 * goes through the center and a point.
42 * So there's such a cone for each position, and a
43 * point is in the solid angle if it lies within any of these cones.
44 * The width of the cones can be adjusted.
46 * The method is implemented by partitioning the surface of the unit sphere
47 * into bins using the polar coordinates \f$(\theta, \phi)\f$.
48 * The partitioning is always uniform in the zenith angle \f$\theta\f$,
49 * while the partitioning in the azimuthal angle \f$\phi\f$ varies.
50 * For each reference point, the unit vector from the center to the point
51 * is constructed, and it is stored in all the bins that overlap with the
52 * cone defined by the point.
53 * Bins that are completely covered by a single cone are marked as such.
54 * Checking whether a point is in the solid angle is then straightforward
55 * with this data structure: one finds the bin that corresponds to the point,
56 * and checks whether the bin is completely covered. If it is not, one
57 * additionally needs to check whether it is within the specified cutoff of
58 * any of the stored points.
60 * The above construction gives quite a lot of flexibility for constructing
61 * the bins without modifying the rest of the code.
62 * The current (quite inefficient) implementation is discussed below, but
63 * it should be optimized to get the most out of the code.
65 * The current way of constructing the bins constructs the boundaries
66 * statically: the bin size in the zenith direction is set to approximately
67 * half the angle cutoff, and the bins in the azimuthal direction have
68 * sizes such that the shortest edge of the bin is approximately equal to
69 * half the angle cutoff (for the regions close to the poles, a single bin
71 * Each reference point is then added to the bins as follows:
72 * -# Find the zenith angle range that is spanned by the cone centered at the
73 * point (this is simple addition/subtraction).
74 * -# Calculate the maximal span of the cone in the azimuthal direction using
76 * \f[\sin \Delta \phi_{max} = \frac{\sin \alpha}{\sin \theta}\f]
77 * (a sine formula in spherical coordinates),
78 * where \f$\alpha\f$ is the width of the cone and \f$\theta\f$ is the
79 * zenith angle of the cone center.
80 * Similarly, the zenith angle at which this extent is achieved is
82 * \f[\cos \theta_{max} = \frac{\cos \theta}{\cos \alpha}\f]
83 * (Pythagoras's theorem in spherical coordinates).
84 * -# For each zenith angle bin that is at least partially covered by the
85 * cone, calculate the span of the cone at the edges using
86 * \f[\sin^2 \frac{\Delta \phi}{2} = \frac{\sin^2 \frac{\alpha}{2} - \sin^2 \frac{\theta - \theta'}{2}}{\sin \theta \sin \theta'}\f]
87 * (distance in spherical geometry),
88 * where \f$\theta'\f$ is the zenith angle of the bin edge.
89 * Treat zenith angle bins that are completely covered by the cone (in the
90 * case that the cone is centered close to the pole) as a special case.
91 * -# Using the values calculated above, loop through the azimuthal bins that
92 * are partially or completely covered by the cone and update them.
94 * The total solid angle (for covered fraction calculations) is estimated by
95 * taking the total area of completely covered bins plus
96 * half the area of partially covered bins.
97 * The second one is an approximation, but should give reasonable estimates
98 * for the averages as well as in cases where the bin size is small.
102 * Implements the \ref sm_insolidangle "insolidangle" selection method.
105 * The implementation could be optimized quite a bit.
108 * Move the covered fraction stuff somewhere else and make it more generic
109 * (along the lines it is handled in selection.h and trajana.h in the old C
112 * \author Teemu Murtola <teemu.murtola@gmail.com>
113 * \ingroup module_selection
121 #include "gromacs/legacyheaders/macros.h"
122 #include "gromacs/math/units.h"
123 #include "gromacs/math/utilities.h"
124 #include "gromacs/math/vec.h"
125 #include "gromacs/pbcutil/pbc.h"
126 #include "gromacs/selection/indexutil.h"
127 #include "gromacs/selection/position.h"
128 #include "gromacs/selection/selection.h"
129 #include "gromacs/utility/exceptions.h"
130 #include "gromacs/utility/smalloc.h"
133 #include "selmethod.h"
140 * Internal data structure for the \p insolidangle selection method.
142 * \see \c t_partition
144 * \ingroup module_selection
148 /** Left edge of the partition. */
150 /** Bin index corresponding to this partition. */
156 * Internal data structure for the \p insolidangle selection method.
158 * Describes the surface partitioning within one slice along the zenith angle.
159 * The slice from azimuthal angle \p p[i].left to \p p[i+1].left belongs to
162 * \ingroup module_selection
166 /** Number of partition items (\p p contains \p n+1 items). */
168 /** Array of partition edges and corresponding bins. */
174 * Internal data structure for the \p insolidangle selection method.
176 * Contains the reference points that partially cover a certain region on the
177 * surface of the unit sphere.
178 * If \p n is -1, the whole region described by the bin is covered.
180 * \ingroup module_selection
184 /** Number of points in the array \p x, -1 if whole bin covered. */
186 /** Number of elements allocated for \p x. */
188 /** Array of points that partially cover the bin. */
190 } t_spheresurfacebin
;
194 * Data structure for the \p insolidangle selection method.
196 * All angle values are in the units of radians.
198 * \ingroup module_selection
202 /** Center of the solid angle. */
203 gmx_ana_pos_t center
;
204 /** Positions that span the solid angle. */
208 /** Estimate of the covered fraction. */
211 /** Cutoff for the cosine (equals cos(angcut)). */
213 /** Bin size to be used as the target bin size when constructing the bins. */
216 /** Number of bins in the \p tbin array. */
218 /** Size of one bin in the zenith angle direction. */
220 /** Array of zenith angle slices. */
222 /** Number of elements allocated for the \p bin array. */
224 /** Number of elements used in the \p bin array. */
226 /** Array of individual bins. */
227 t_spheresurfacebin
*bin
;
228 } t_methoddata_insolidangle
;
231 * Allocates data for the \p insolidangle selection method.
233 * \param[in] npar Not used (should be 3).
234 * \param[in,out] param Method parameters (should point to
235 * \ref smparams_insolidangle).
236 * \returns Pointer to the allocated data (\ref t_methoddata_insolidangle).
238 * Allocates memory for a \ref t_methoddata_insolidangle structure and
239 * initializes the parameter as follows:
240 * - \p center defines the value for t_methoddata_insolidangle::center.
241 * - \p span defines the value for t_methoddata_insolidangle::span.
242 * - \p cutoff defines the value for t_methoddata_insolidangle::angcut.
245 init_data_insolidangle(int npar
, gmx_ana_selparam_t
*param
);
247 * Initializes the \p insolidangle selection method.
249 * \param top Not used.
250 * \param npar Not used.
251 * \param param Not used.
252 * \param data Pointer to \ref t_methoddata_insolidangle to initialize.
253 * \returns 0 on success, -1 on failure.
255 * Converts t_methoddata_insolidangle::angcut to radians and allocates
256 * and allocates memory for the bins used during the evaluation.
259 init_insolidangle(t_topology
* top
, int npar
, gmx_ana_selparam_t
* param
, void *data
);
260 /** Frees the data allocated for the \p insolidangle selection method. */
262 free_data_insolidangle(void *data
);
264 * Initializes the evaluation of the \p insolidangle selection method for a frame.
266 * \param[in] top Not used.
267 * \param[in] fr Not used.
268 * \param[in] pbc PBC structure.
269 * \param data Should point to a \ref t_methoddata_insolidangle.
271 * Creates a lookup structure that enables fast queries of whether a point
272 * is within the solid angle or not.
275 init_frame_insolidangle(t_topology
* top
, t_trxframe
* fr
, t_pbc
*pbc
, void *data
);
276 /** Internal helper function for evaluate_insolidangle(). */
278 accept_insolidangle(rvec x
, t_pbc
*pbc
, void *data
);
279 /** Evaluates the \p insolidangle selection method. */
281 evaluate_insolidangle(t_topology
* /* top */, t_trxframe
* /* fr */, t_pbc
*pbc
,
282 gmx_ana_pos_t
*pos
, gmx_ana_selvalue_t
*out
, void *data
);
284 /** Calculates the distance between unit vectors. */
286 sph_distc(rvec x1
, rvec x2
);
287 /** Does a binary search on a \p t_partition to find a bin for a value. */
289 find_partition_bin(t_partition
*p
, real value
);
290 /** Finds a bin that corresponds to a location on the unit sphere surface. */
292 find_surface_bin(t_methoddata_insolidangle
*surf
, rvec x
);
293 /** Clears/initializes the bins on the unit sphere surface. */
295 clear_surface_points(t_methoddata_insolidangle
*surf
);
296 /** Frees memory allocated for storing the reference points in the surface bins. */
298 free_surface_points(t_methoddata_insolidangle
*surf
);
299 /** Adds a reference point to a given bin. */
301 add_surface_point(t_methoddata_insolidangle
*surf
, int tbin
, int pbin
, rvec x
);
302 /** Marks a bin as completely covered. */
304 mark_surface_covered(t_methoddata_insolidangle
*surf
, int tbin
, int pbin
);
305 /** Helper function for store_surface_point() to update a single zenith angle bin. */
307 update_surface_bin(t_methoddata_insolidangle
*surf
, int tbin
,
308 real phi
, real pdelta1
, real pdelta2
, real pdeltamax
,
310 /** Adds a single reference point and updates the surface bins. */
312 store_surface_point(t_methoddata_insolidangle
*surf
, rvec x
);
314 * Optimizes the surface bins for faster searching.
316 * \param[in,out] surf Surface data structure.
318 * Currently, this function does nothing.
321 optimize_surface_points(t_methoddata_insolidangle
*surf
);
322 /** Estimates the area covered by the reference cones. */
324 estimate_covered_fraction(t_methoddata_insolidangle
*surf
);
325 /** Checks whether a point lies within a solid angle. */
327 is_surface_covered(t_methoddata_insolidangle
*surf
, rvec x
);
329 /** Parameters for the \p insolidangle selection method. */
330 static gmx_ana_selparam_t smparams_insolidangle
[] = {
331 {"center", {POS_VALUE
, 1, {NULL
}}, NULL
, SPAR_DYNAMIC
},
332 {"span", {POS_VALUE
, -1, {NULL
}}, NULL
, SPAR_DYNAMIC
| SPAR_VARNUM
},
333 {"cutoff", {REAL_VALUE
, 1, {NULL
}}, NULL
, SPAR_OPTIONAL
},
336 /** Help text for the \p insolidangle selection method. */
337 static const char *const help_insolidangle
[] = {
340 " insolidangle center POS span POS_EXPR [cutoff REAL]",
342 "This keyword selects atoms that are within [TT]REAL[tt] degrees",
343 "(default=5) of any position in [TT]POS_EXPR[tt] as seen from [TT]POS[tt]",
344 "a position expression that evaluates to a single position), i.e., atoms",
345 "in the solid angle spanned by the positions in [TT]POS_EXPR[tt] and",
346 "centered at [TT]POS[tt].[PAR]"
348 "Technically, the solid angle is constructed as a union of small cones",
349 "whose tip is at [TT]POS[tt] and the axis goes through a point in",
350 "[TT]POS_EXPR[tt]. There is such a cone for each position in",
351 "[TT]POS_EXPR[tt], and point is in the solid angle if it lies within any",
352 "of these cones. The cutoff determines the width of the cones.",
355 /** Selection method data for the \p insolidangle method. */
356 gmx_ana_selmethod_t sm_insolidangle
= {
357 "insolidangle", GROUP_VALUE
, SMETH_DYNAMIC
,
358 asize(smparams_insolidangle
), smparams_insolidangle
,
359 &init_data_insolidangle
,
363 &free_data_insolidangle
,
364 &init_frame_insolidangle
,
366 &evaluate_insolidangle
,
367 {"insolidangle center POS span POS_EXPR [cutoff REAL]",
368 "Selecting atoms in a solid angle",
369 asize(help_insolidangle
), help_insolidangle
},
373 init_data_insolidangle(int /* npar */, gmx_ana_selparam_t
*param
)
375 t_methoddata_insolidangle
*data
= new t_methoddata_insolidangle();
379 data
->distccut
= 0.0;
380 data
->targetbinsize
= 0.0;
383 data
->tbinsize
= 0.0;
389 param
[0].val
.u
.p
= &data
->center
;
390 param
[1].val
.u
.p
= &data
->span
;
391 param
[2].val
.u
.r
= &data
->angcut
;
396 init_insolidangle(t_topology
* /* top */, int /* npar */, gmx_ana_selparam_t
* /* param */, void *data
)
398 t_methoddata_insolidangle
*surf
= (t_methoddata_insolidangle
*)data
;
401 if (surf
->angcut
<= 0)
403 GMX_THROW(gmx::InvalidInputError("Angle cutoff should be > 0"));
406 surf
->angcut
*= DEG2RAD
;
408 surf
->distccut
= -cos(surf
->angcut
);
409 surf
->targetbinsize
= surf
->angcut
/ 2;
410 surf
->ntbins
= static_cast<int>(M_PI
/ surf
->targetbinsize
);
411 surf
->tbinsize
= (180.0 / surf
->ntbins
)*DEG2RAD
;
413 snew(surf
->tbin
, static_cast<int>(M_PI
/ surf
->tbinsize
) + 1);
415 for (i
= 0; i
< surf
->ntbins
; ++i
)
417 c
= static_cast<int>(max(sin(surf
->tbinsize
*i
),
418 sin(surf
->tbinsize
*(i
+1)))
419 * M_2PI
/ surf
->targetbinsize
) + 1;
420 snew(surf
->tbin
[i
].p
, c
+1);
424 snew(surf
->bin
, surf
->maxbins
);
428 * \param data Data to free (should point to a \ref t_methoddata_insolidangle).
430 * Frees the memory allocated for \c t_methoddata_insolidangle::center and
431 * \c t_methoddata_insolidangle::span, as well as the memory for the internal
435 free_data_insolidangle(void *data
)
437 t_methoddata_insolidangle
*d
= (t_methoddata_insolidangle
*)data
;
442 for (i
= 0; i
< d
->ntbins
; ++i
)
448 free_surface_points(d
);
454 init_frame_insolidangle(t_topology
* /* top */, t_trxframe
* /* fr */, t_pbc
*pbc
, void *data
)
456 t_methoddata_insolidangle
*d
= (t_methoddata_insolidangle
*)data
;
460 free_surface_points(d
);
461 clear_surface_points(d
);
462 for (i
= 0; i
< d
->span
.count(); ++i
)
466 pbc_dx(pbc
, d
->span
.x
[i
], d
->center
.x
[0], dx
);
470 rvec_sub(d
->span
.x
[i
], d
->center
.x
[0], dx
);
473 store_surface_point(d
, dx
);
475 optimize_surface_points(d
);
480 * \param[in] x Test point.
481 * \param[in] pbc PBC data (if NULL, no PBC are used).
482 * \param[in] data Pointer to a \c t_methoddata_insolidangle data structure.
483 * \returns true if \p x is within the solid angle, false otherwise.
486 accept_insolidangle(rvec x
, t_pbc
*pbc
, void *data
)
488 t_methoddata_insolidangle
*d
= (t_methoddata_insolidangle
*)data
;
493 pbc_dx(pbc
, x
, d
->center
.x
[0], dx
);
497 rvec_sub(x
, d
->center
.x
[0], dx
);
500 return is_surface_covered(d
, dx
);
504 * See sel_updatefunc() for description of the parameters.
505 * \p data should point to a \c t_methoddata_insolidangle.
507 * Calculates which atoms in \p g are within the solid angle spanned by
508 * \c t_methoddata_insolidangle::span and centered at
509 * \c t_methoddata_insolidangle::center, and stores the result in \p out->u.g.
512 evaluate_insolidangle(t_topology
* /* top */, t_trxframe
* /* fr */, t_pbc
*pbc
,
513 gmx_ana_pos_t
*pos
, gmx_ana_selvalue_t
*out
, void *data
)
516 for (int b
= 0; b
< pos
->count(); ++b
)
518 if (accept_insolidangle(pos
->x
[b
], pbc
, data
))
520 gmx_ana_pos_add_to_group(out
->u
.g
, pos
, b
);
526 * \param[in] sel Selection element to query.
527 * \returns true if the covered fraction can be estimated for \p sel with
528 * _gmx_selelem_estimate_coverfrac(), false otherwise.
531 _gmx_selelem_can_estimate_cover(const gmx::SelectionTreeElement
&sel
)
533 if (sel
.type
== SEL_BOOLEAN
&& sel
.u
.boolt
== BOOL_OR
)
538 bool bDynFound
= false;
539 gmx::SelectionTreeElementPointer child
= sel
.child
;
542 if (child
->type
== SEL_EXPRESSION
)
544 if (child
->u
.expr
.method
->name
== sm_insolidangle
.name
)
546 if (bFound
|| bDynFound
)
552 else if (child
->u
.expr
.method
553 && (child
->u
.expr
.method
->flags
& SMETH_DYNAMIC
))
562 else if (!_gmx_selelem_can_estimate_cover(*child
))
572 * \param[in] sel Selection for which the fraction should be calculated.
573 * \returns Fraction of angles covered by the selection (between zero and one).
575 * The return value is undefined if _gmx_selelem_can_estimate_cover() returns
577 * Should be called after gmx_ana_evaluate_selections() has been called for the
581 _gmx_selelem_estimate_coverfrac(const gmx::SelectionTreeElement
&sel
)
585 if (sel
.type
== SEL_EXPRESSION
&& sel
.u
.expr
.method
->name
== sm_insolidangle
.name
)
587 t_methoddata_insolidangle
*d
= (t_methoddata_insolidangle
*)sel
.u
.expr
.mdata
;
590 d
->cfrac
= estimate_covered_fraction(d
);
594 if (sel
.type
== SEL_BOOLEAN
&& sel
.u
.boolt
== BOOL_NOT
)
596 cfrac
= _gmx_selelem_estimate_coverfrac(*sel
.child
);
604 /* Here, we assume that the selection is simple enough */
605 gmx::SelectionTreeElementPointer child
= sel
.child
;
608 cfrac
= _gmx_selelem_estimate_coverfrac(*child
);
619 * \param[in] x1 Unit vector 1.
620 * \param[in] x2 Unit vector 2.
621 * \returns Minus the dot product of \p x1 and \p x2.
623 * This function is used internally to calculate the distance between the
624 * unit vectors \p x1 and \p x2 to find out whether \p x2 is within the
625 * cone centered at \p x1. Currently, the cosine of the angle is used
626 * for efficiency, and the minus is there to make it behave like a normal
627 * distance (larger values mean longer distances).
630 sph_distc(rvec x1
, rvec x2
)
632 return -iprod(x1
, x2
);
636 * \param[in] p Partition to search.
637 * \param[in] value Value to search for.
638 * \returns The partition index in \p p that contains \p value.
640 * If \p value is outside the range of \p p, the first/last index is returned.
641 * Otherwise, the return value \c i satisfies \c p->p[i].left<=value and
642 * \c p->p[i+1].left>value
645 find_partition_bin(t_partition
*p
, real value
)
647 int pmin
, pmax
, pbin
;
649 /* Binary search the partition */
650 pmin
= 0; pmax
= p
->n
;
651 while (pmax
> pmin
+ 1)
653 pbin
= pmin
+ (pmax
- pmin
) / 2;
654 if (p
->p
[pbin
].left
<= value
)
668 * \param[in] surf Surface data structure to search.
669 * \param[in] x Unit vector to find.
670 * \returns The bin index that contains \p x.
672 * The return value is an index to the \p surf->bin array.
675 find_surface_bin(t_methoddata_insolidangle
*surf
, rvec x
)
681 phi
= atan2(x
[YY
], x
[XX
]);
682 tbin
= static_cast<int>(floor(theta
/ surf
->tbinsize
));
683 if (tbin
>= surf
->ntbins
)
685 tbin
= surf
->ntbins
- 1;
687 pbin
= find_partition_bin(&surf
->tbin
[tbin
], phi
);
688 return surf
->tbin
[tbin
].p
[pbin
].bin
;
692 * \param[in,out] surf Surface data structure.
694 * Clears the reference points from the bins and (re)initializes the edges
695 * of the azimuthal bins.
698 clear_surface_points(t_methoddata_insolidangle
*surf
)
703 for (i
= 0; i
< surf
->ntbins
; ++i
)
705 c
= static_cast<int>(min(sin(surf
->tbinsize
*i
),
706 sin(surf
->tbinsize
*(i
+1)))
707 * M_2PI
/ surf
->targetbinsize
) + 1;
713 for (j
= 0; j
< c
; ++j
)
715 surf
->tbin
[i
].p
[j
].left
= -M_PI
+ j
*M_2PI
/c
- 0.0001;
716 surf
->tbin
[i
].p
[j
].bin
= surf
->nbins
;
717 surf
->bin
[surf
->nbins
].n
= 0;
720 surf
->tbin
[i
].p
[c
].left
= M_PI
+ 0.0001;
721 surf
->tbin
[i
].p
[c
].bin
= -1;
726 * \param[in,out] surf Surface data structure.
729 free_surface_points(t_methoddata_insolidangle
*surf
)
733 for (i
= 0; i
< surf
->nbins
; ++i
)
737 sfree(surf
->bin
[i
].x
);
739 surf
->bin
[i
].n_alloc
= 0;
740 surf
->bin
[i
].x
= NULL
;
745 * \param[in,out] surf Surface data structure.
746 * \param[in] tbin Bin number in the zenith angle direction.
747 * \param[in] pbin Bin number in the azimuthal angle direction.
748 * \param[in] x Point to store.
751 add_surface_point(t_methoddata_insolidangle
*surf
, int tbin
, int pbin
, rvec x
)
755 bin
= surf
->tbin
[tbin
].p
[pbin
].bin
;
756 /* Return if bin is already completely covered */
757 if (surf
->bin
[bin
].n
== -1)
761 /* Allocate more space if necessary */
762 if (surf
->bin
[bin
].n
== surf
->bin
[bin
].n_alloc
)
764 surf
->bin
[bin
].n_alloc
+= 10;
765 srenew(surf
->bin
[bin
].x
, surf
->bin
[bin
].n_alloc
);
767 /* Add the point to the bin */
768 copy_rvec(x
, surf
->bin
[bin
].x
[surf
->bin
[bin
].n
]);
773 * \param[in,out] surf Surface data structure.
774 * \param[in] tbin Bin number in the zenith angle direction.
775 * \param[in] pbin Bin number in the azimuthal angle direction.
778 mark_surface_covered(t_methoddata_insolidangle
*surf
, int tbin
, int pbin
)
782 bin
= surf
->tbin
[tbin
].p
[pbin
].bin
;
783 surf
->bin
[bin
].n
= -1;
787 * \param[in,out] surf Surface data structure.
788 * \param[in] tbin Bin number in the zenith angle direction.
789 * \param[in] phi Azimuthal angle of \p x.
790 * \param[in] pdelta1 Width of the cone at the lower edge of \p tbin.
791 * \param[in] pdelta2 Width of the cone at the uppper edge of \p tbin.
792 * \param[in] pdeltamax Max. width of the cone inside \p tbin.
793 * \param[in] x Point to store (should have unit length).
796 update_surface_bin(t_methoddata_insolidangle
*surf
, int tbin
,
797 real phi
, real pdelta1
, real pdelta2
, real pdeltamax
,
800 real pdelta
, phi1
, phi2
;
801 int pbin1
, pbin2
, pbiniter
, pbin
;
803 /* Find the edges of the bins affected */
804 pdelta
= max(max(pdelta1
, pdelta2
), pdeltamax
);
808 pbin
= find_partition_bin(&surf
->tbin
[tbin
], phi1
);
813 pbin
= find_partition_bin(&surf
->tbin
[tbin
], phi1
+ M_2PI
);
814 pbin1
= pbin
- surf
->tbin
[tbin
].n
;
819 pbin2
= find_partition_bin(&surf
->tbin
[tbin
], phi2
);
823 pbin2
= find_partition_bin(&surf
->tbin
[tbin
], phi2
- M_2PI
);
824 pbin2
+= surf
->tbin
[tbin
].n
;
827 if (pbin2
- pbin1
> surf
->tbin
[tbin
].n
)
829 pbin2
= pbin1
+ surf
->tbin
[tbin
].n
;
831 /* Find the edges of completely covered region */
832 pdelta
= min(pdelta1
, pdelta2
);
839 /* Loop over all affected bins */
840 for (pbiniter
= pbin1
; pbiniter
!= pbin2
; ++pbiniter
, ++pbin
)
842 /* Wrap bin around if end reached */
843 if (pbin
== surf
->tbin
[tbin
].n
)
849 /* Check if bin is completely covered and update */
850 if (surf
->tbin
[tbin
].p
[pbin
].left
>= phi1
851 && surf
->tbin
[tbin
].p
[pbin
+1].left
<= phi2
)
853 mark_surface_covered(surf
, tbin
, pbin
);
857 add_surface_point(surf
, tbin
, pbin
, x
);
863 * \param[in,out] surf Surface data structure.
864 * \param[in] x Point to store (should have unit length).
866 * Finds all the bins covered by the cone centered at \p x and calls
867 * update_surface_bin() to update them.
870 store_surface_point(t_methoddata_insolidangle
*surf
, rvec x
)
873 real pdeltamax
, tmax
;
874 real theta1
, theta2
, pdelta1
, pdelta2
;
878 phi
= atan2(x
[YY
], x
[XX
]);
879 /* Find the maximum extent in the phi direction */
880 if (theta
<= surf
->angcut
)
885 else if (theta
>= M_PI
- surf
->angcut
)
892 pdeltamax
= asin(sin(surf
->angcut
) / sin(theta
));
893 tmax
= acos(cos(theta
) / cos(surf
->angcut
));
895 /* Find the first affected bin */
896 tbin
= max(static_cast<int>(floor((theta
- surf
->angcut
) / surf
->tbinsize
)), 0);
897 theta1
= tbin
* surf
->tbinsize
;
898 if (theta1
< theta
- surf
->angcut
)
906 /* Loop through all affected bins */
907 while (tbin
< ceil((theta
+ surf
->angcut
) / surf
->tbinsize
)
908 && tbin
< surf
->ntbins
)
910 /* Calculate the next boundaries */
911 theta2
= (tbin
+1) * surf
->tbinsize
;
912 if (theta2
> theta
+ surf
->angcut
)
914 /* The circle is completely outside the cone */
917 else if (theta2
<= -(theta
- surf
->angcut
)
918 || theta2
>= M_2PI
- (theta
+ surf
->angcut
)
919 || tbin
== surf
->ntbins
- 1)
921 /* The circle is completely inside the cone, or we are in the
922 * 360 degree bin covering the pole. */
927 /* TODO: This formula is numerically unstable if theta is very
928 * close to the pole. In practice, it probably does not matter
929 * much, but it would be nicer to adjust the theta bin boundaries
930 * such that the case above catches this instead of falling through
932 pdelta2
= 2*asin(sqrt(
933 (sqr(sin(surf
->angcut
/2)) - sqr(sin((theta2
-theta
)/2))) /
934 (sin(theta
) * sin(theta2
))));
937 if (tmax
>= theta1
&& tmax
<= theta2
)
939 update_surface_bin(surf
, tbin
, phi
, pdelta1
, pdelta2
, pdeltamax
, x
);
943 update_surface_bin(surf
, tbin
, phi
, pdelta1
, pdelta2
, 0, x
);
953 optimize_surface_points(t_methoddata_insolidangle
* /* surf */)
955 /* TODO: Implement */
959 * \param[in] surf Surface data structure.
960 * \returns An estimate for the area covered by the reference points.
963 estimate_covered_fraction(t_methoddata_insolidangle
*surf
)
966 real cfrac
, tfrac
, pfrac
;
969 for (t
= 0; t
< surf
->ntbins
; ++t
)
971 tfrac
= cos(t
* surf
->tbinsize
) - cos((t
+1) * surf
->tbinsize
);
972 for (p
= 0; p
< surf
->tbin
[t
].n
; ++p
)
974 pfrac
= surf
->tbin
[t
].p
[p
+1].left
- surf
->tbin
[t
].p
[p
].left
;
975 n
= surf
->bin
[surf
->tbin
[t
].p
[p
].bin
].n
;
976 if (n
== -1) /* Bin completely covered */
978 cfrac
+= tfrac
* pfrac
;
980 else if (n
> 0) /* Bin partially covered */
982 cfrac
+= tfrac
* pfrac
/ 2; /* A rough estimate */
986 return cfrac
/ (4*M_PI
);
990 * \param[in] surf Surface data structure to search.
991 * \param[in] x Unit vector to check.
992 * \returns true if \p x is within the solid angle, false otherwise.
995 is_surface_covered(t_methoddata_insolidangle
*surf
, rvec x
)
999 bin
= find_surface_bin(surf
, x
);
1000 /* Check for completely covered bin */
1001 if (surf
->bin
[bin
].n
== -1)
1005 /* Check each point that partially covers the bin */
1006 for (i
= 0; i
< surf
->bin
[bin
].n
; ++i
)
1008 if (sph_distc(x
, surf
->bin
[bin
].x
[i
]) < surf
->distccut
)