3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
9 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
10 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
11 * Copyright (c) 2001-2009, The GROMACS development team,
12 * check out http://www.gromacs.org for more information.
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * as published by the Free Software Foundation; either version 2
17 * of the License, or (at your option) any later version.
19 * If you want to redistribute modifications, please consider that
20 * scientific software is very special. Version control is crucial -
21 * bugs must be traceable. We will be happy to consider code for
22 * inclusion in the official distribution, but derived work must not
23 * be called official GROMACS. Details are found in the README & COPYING
24 * files - if they are missing, get the official version at www.gromacs.org.
26 * To help us fund GROMACS development, we humbly ask that you cite
27 * the papers on the package - you can find them in the top README file.
29 * For more info, check our website at http://www.gromacs.org
31 /*! \page page_module_selection_insolidangle Selection method: insolidangle
33 * This method selects a subset of particles that are located in a solid
34 * angle defined by a center and a set of points.
35 * The solid angle is constructed as a union of small cones whose axis
36 * goes through the center and a point.
37 * So there's such a cone for each position, and a
38 * point is in the solid angle if it lies within any of these cones.
39 * The width of the cones can be adjusted.
43 * The method is implemented by partitioning the surface of the unit sphere
44 * into bins using the polar coordinates \f$(\theta, \phi)\f$.
45 * The partitioning is always uniform in the zenith angle \f$\theta\f$,
46 * while the partitioning in the azimuthal angle \f$\phi\f$ varies.
47 * For each reference point, the unit vector from the center to the point
48 * is constructed, and it is stored in all the bins that overlap with the
49 * cone defined by the point.
50 * Bins that are completely covered by a single cone are marked as such.
51 * Checking whether a point is in the solid angle is then straightforward
52 * with this data structure: one finds the bin that corresponds to the point,
53 * and checks whether the bin is completely covered. If it is not, one
54 * additionally needs to check whether it is within the specified cutoff of
55 * any of the stored points.
57 * The above construction gives quite a lot of flexibility for constructing
58 * the bins without modifying the rest of the code.
59 * The current (quite inefficient) implementation is discussed below, but
60 * it should be optimized to get the most out of the code.
62 * The current way of constructing the bins constructs the boundaries
63 * statically: the bin size in the zenith direction is set to approximately
64 * half the angle cutoff, and the bins in the azimuthal direction have
65 * sizes such that the shortest edge of the bin is approximately equal to
66 * half the angle cutoff (for the regions close to the poles, a single bin
68 * Each reference point is then added to the bins as follows:
69 * -# Find the zenith angle range that is spanned by the cone centered at the
70 * point (this is simple addition/subtraction).
71 * -# Calculate the maximal span of the cone in the azimuthal direction using
73 * \f[\sin \Delta \phi_{max} = \frac{\sin \alpha}{\sin \theta}\f]
74 * (a sine formula in spherical coordinates),
75 * where \f$\alpha\f$ is the width of the cone and \f$\theta\f$ is the
76 * zenith angle of the cone center.
77 * Similarly, the zenith angle at which this extent is achieved is
79 * \f[\cos \theta_{max} = \frac{\cos \theta}{\cos \alpha}\f]
80 * (Pythagoras's theorem in spherical coordinates).
81 * -# For each zenith angle bin that is at least partially covered by the
82 * cone, calculate the span of the cone at the edges using
83 * \f[\sin^2 \frac{\Delta \phi}{2} = \frac{\sin^2 \frac{\alpha}{2} - \sin^2 \frac{\theta - \theta'}{2}}{\sin \theta \sin \theta'}\f]
84 * (distance in spherical geometry),
85 * where \f$\theta'\f$ is the zenith angle of the bin edge.
86 * -# Using the values calculated above, loop through the azimuthal bins that
87 * are partially or completely covered by the cone and update them.
89 * The total solid angle (for covered fraction calculations) is estimated by
90 * taking the total area of completely covered bins plus
91 * half the area of partially covered bins.
92 * The second one is an approximation, but should give reasonable estimates
93 * for the averages as well as in cases where the bin size is small.
97 * Implements the \ref sm_insolidangle "insolidangle" selection method.
100 * The implementation could be optimized quite a bit.
103 * Move the covered fraction stuff somewhere else and make it more generic
104 * (along the lines it is handled in selection.h and trajana.h in the old C
107 * \author Teemu Murtola <teemu.murtola@cbr.su.se>
108 * \ingroup module_selection
123 // FIXME: Should really be in the beginning, but causes compilation errors
126 #include "gromacs/fatalerror/exceptions.h"
127 #include "gromacs/selection/indexutil.h"
128 #include "gromacs/selection/position.h"
129 #include "gromacs/selection/selection.h"
130 #include "gromacs/selection/selmethod.h"
138 * Internal data structure for the \p insolidangle selection method.
140 * \see \c t_partition
144 /** Left edge of the partition. */
146 /** Bin index corresponding to this partition. */
151 * Internal data structure for the \p insolidangle selection method.
153 * Describes the surface partitioning within one slice along the zenith angle.
154 * The slice from azimuthal angle \p p[i].left to \p p[i+1].left belongs to
159 /** Number of partition items (\p p contains \p n+1 items). */
161 /** Array of partition edges and corresponding bins. */
166 * Internal data structure for the \p insolidangle selection method.
168 * Contains the reference points that partially cover a certain region on the
169 * surface of the unit sphere.
170 * If \p n is -1, the whole region described by the bin is covered.
174 /** Number of points in the array \p x, -1 if whole bin covered. */
176 /** Number of elements allocated for \p x. */
178 /** Array of points that partially cover the bin. */
180 } t_spheresurfacebin;
183 * Data structure for the \p insolidangle selection method.
185 * All angle values are in the units of radians.
189 /** Center of the solid angle. */
190 gmx_ana_pos_t center;
191 /** Positions that span the solid angle. */
195 /** Estimate of the covered fraction. */
198 /** Cutoff for the cosine (equals cos(angcut)). */
200 /** Bin size to be used as the target bin size when constructing the bins. */
203 /** Number of bins in the \p tbin array. */
205 /** Size of one bin in the zenith angle direction. */
207 /** Array of zenith angle slices. */
209 /** Number of elements allocated for the \p bin array. */
211 /** Number of elements used in the \p bin array. */
213 /** Array of individual bins. */
214 t_spheresurfacebin *bin;
215 } t_methoddata_insolidangle;
217 /** Allocates data for the \p insolidangle selection method. */
219 init_data_insolidangle(int npar, gmx_ana_selparam_t *param);
220 /** Initializes the \p insolidangle selection method. */
222 init_insolidangle(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
223 /** Sets the COM/COG data for the \p insolidangle selection method. */
225 set_comg_insolidangle(gmx_ana_pos_t *pos, void *data);
226 /** Frees the data allocated for the \p insolidangle selection method. */
228 free_data_insolidangle(void *data);
229 /** Initializes the evaluation of the \p insolidangle selection method for a frame. */
231 init_frame_insolidangle(t_topology *top, t_trxframe *fr, t_pbc *pbc, void *data);
232 /** Internal helper function for evaluate_insolidangle(). */
234 accept_insolidangle(rvec x, t_pbc *pbc, void *data);
235 /** Evaluates the \p insolidangle selection method. */
237 evaluate_insolidangle(t_topology *top, t_trxframe *fr, t_pbc *pbc,
238 gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void *data);
240 /** Calculates the distance between unit vectors. */
242 sph_distc(rvec x1, rvec x2);
243 /** Does a binary search on a \p t_partition to find a bin for a value. */
245 find_partition_bin(t_partition *p, real value);
246 /** Finds a bin that corresponds to a location on the unit sphere surface. */
248 find_surface_bin(t_methoddata_insolidangle *surf, rvec x);
249 /** Clears/initializes the bins on the unit sphere surface. */
251 clear_surface_points(t_methoddata_insolidangle *surf);
252 /** Frees memory allocated for storing the reference points in the surface bins. */
254 free_surface_points(t_methoddata_insolidangle *surf);
255 /** Adds a reference point to a given bin. */
257 add_surface_point(t_methoddata_insolidangle *surf, int tbin, int pbin, rvec x);
258 /** Marks a bin as completely covered. */
260 mark_surface_covered(t_methoddata_insolidangle *surf, int tbin, int pbin);
261 /** Helper function for store_surface_point() to update a single zenith angle bin. */
263 update_surface_bin(t_methoddata_insolidangle *surf, int tbin,
264 real phi, real pdelta1, real pdelta2, real pdeltamax,
266 /** Adds a single reference point and updates the surface bins. */
268 store_surface_point(t_methoddata_insolidangle *surf, rvec x);
269 /** Optimizes the surface bins for faster searching. */
271 optimize_surface_points(t_methoddata_insolidangle *surf);
272 /** Estimates the area covered by the reference cones. */
274 estimate_covered_fraction(t_methoddata_insolidangle *surf);
275 /** Checks whether a point lies within a solid angle. */
277 is_surface_covered(t_methoddata_insolidangle *surf, rvec x);
279 /** Parameters for the \p insolidangle selection method. */
280 static gmx_ana_selparam_t smparams_insolidangle[] = {
281 {"center", {POS_VALUE, 1, {NULL}}, NULL, SPAR_DYNAMIC},
282 {"span", {POS_VALUE, -1, {NULL}}, NULL, SPAR_DYNAMIC | SPAR_VARNUM},
283 {"cutoff", {REAL_VALUE, 1, {NULL}}, NULL, SPAR_OPTIONAL},
286 /** Help text for the \p insolidangle selection method. */
287 static const char *help_insolidangle[] = {
288 "SELECTING ATOMS IN A SOLID ANGLE[PAR]",
290 "[TT]insolidangle center POS span POS_EXPR [cutoff REAL][tt][PAR]",
292 "This keyword selects atoms that are within [TT]REAL[tt] degrees",
293 "(default=5) of any position in [TT]POS_EXPR[tt] as seen from [TT]POS[tt]",
294 "a position expression that evaluates to a single position), i.e., atoms",
295 "in the solid angle spanned by the positions in [TT]POS_EXPR[tt] and",
296 "centered at [TT]POS[tt].[PAR]"
298 "Technically, the solid angle is constructed as a union of small cones",
299 "whose tip is at [TT]POS[tt] and the axis goes through a point in",
300 "[TT]POS_EXPR[tt]. There is such a cone for each position in",
301 "[TT]POS_EXPR[tt], and point is in the solid angle if it lies within any",
302 "of these cones. The cutoff determines the width of the cones.",
305 /** \internal Selection method data for the \p insolidangle method. */
306 gmx_ana_selmethod_t sm_insolidangle = {
307 "insolidangle", GROUP_VALUE, SMETH_DYNAMIC,
308 asize(smparams_insolidangle), smparams_insolidangle,
309 &init_data_insolidangle,
313 &free_data_insolidangle,
314 &init_frame_insolidangle,
316 &evaluate_insolidangle,
317 {"insolidangle center POS span POS_EXPR [cutoff REAL]",
318 asize(help_insolidangle), help_insolidangle},
322 * \param[in] npar Not used (should be 3).
323 * \param[in,out] param Method parameters (should point to
324 * \ref smparams_insolidangle).
325 * \returns Pointer to the allocated data (\ref t_methoddata_insolidangle).
327 * Allocates memory for a \ref t_methoddata_insolidangle structure and
328 * initializes the parameter as follows:
329 * - \p center defines the value for t_methoddata_insolidangle::center.
330 * - \p span defines the value for t_methoddata_insolidangle::span.
331 * - \p cutoff defines the value for t_methoddata_insolidangle::angcut.
334 init_data_insolidangle(int npar, gmx_ana_selparam_t *param)
336 t_methoddata_insolidangle *data;
340 param[0].val.u.p = &data->center;
341 param[1].val.u.p = &data->span;
342 param[2].val.u.r = &data->angcut;
347 * \param top Not used.
348 * \param npar Not used.
349 * \param param Not used.
350 * \param data Pointer to \ref t_methoddata_insolidangle to initialize.
351 * \returns 0 on success, -1 on failure.
353 * Converts t_methoddata_insolidangle::angcut to radians and allocates
354 * and allocates memory for the bins used during the evaluation.
357 init_insolidangle(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data)
359 t_methoddata_insolidangle *surf = (t_methoddata_insolidangle *)data;
362 if (surf->angcut <= 0)
364 GMX_THROW(gmx::InvalidInputError("Angle cutoff should be > 0"));
367 surf->angcut *= DEG2RAD;
369 surf->distccut = -cos(surf->angcut);
370 surf->targetbinsize = surf->angcut / 2;
371 surf->ntbins = (int) (M_PI / surf->targetbinsize);
372 surf->tbinsize = (180.0 / surf->ntbins)*DEG2RAD;
374 snew(surf->tbin, (int)(M_PI/surf->tbinsize) + 1);
376 for (i = 0; i < surf->ntbins; ++i)
378 c = max(sin(surf->tbinsize*i), sin(surf->tbinsize*(i+1)))
379 * M_2PI / surf->targetbinsize + 1;
380 snew(surf->tbin[i].p, c+1);
384 snew(surf->bin, surf->maxbins);
388 * \param data Data to free (should point to a \ref t_methoddata_insolidangle).
390 * Frees the memory allocated for \c t_methoddata_insolidangle::center and
391 * \c t_methoddata_insolidangle::span, as well as the memory for the internal
395 free_data_insolidangle(void *data)
397 t_methoddata_insolidangle *d = (t_methoddata_insolidangle *)data;
402 for (i = 0; i < d->ntbins; ++i)
408 free_surface_points(d);
413 * \param[in] top Not used.
414 * \param[in] fr Current frame.
415 * \param[in] pbc PBC structure.
416 * \param data Should point to a \ref t_methoddata_insolidangle.
418 * Creates a lookup structure that enables fast queries of whether a point
419 * is within the solid angle or not.
422 init_frame_insolidangle(t_topology *top, t_trxframe *fr, t_pbc *pbc, void *data)
424 t_methoddata_insolidangle *d = (t_methoddata_insolidangle *)data;
428 free_surface_points(d);
429 clear_surface_points(d);
430 for (i = 0; i < d->span.nr; ++i)
434 pbc_dx(pbc, d->span.x[i], d->center.x[0], dx);
438 rvec_sub(d->span.x[i], d->center.x[0], dx);
441 store_surface_point(d, dx);
443 optimize_surface_points(d);
448 * \param[in] x Test point.
449 * \param[in] pbc PBC data (if NULL, no PBC are used).
450 * \param[in] data Pointer to a \c t_methoddata_insolidangle data structure.
451 * \returns TRUE if \p x is within the solid angle, FALSE otherwise.
454 accept_insolidangle(rvec x, t_pbc *pbc, void *data)
456 t_methoddata_insolidangle *d = (t_methoddata_insolidangle *)data;
461 pbc_dx(pbc, x, d->center.x[0], dx);
465 rvec_sub(x, d->center.x[0], dx);
468 return is_surface_covered(d, dx);
472 * See sel_updatefunc() for description of the parameters.
473 * \p data should point to a \c t_methoddata_insolidangle.
475 * Calculates which atoms in \p g are within the solid angle spanned by
476 * \c t_methoddata_insolidangle::span and centered at
477 * \c t_methoddata_insolidangle::center, and stores the result in \p out->u.g.
480 evaluate_insolidangle(t_topology *top, t_trxframe *fr, t_pbc *pbc,
481 gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void *data)
483 t_methoddata_insolidangle *d = (t_methoddata_insolidangle *)data;
487 for (b = 0; b < pos->nr; ++b)
489 if (accept_insolidangle(pos->x[b], pbc, data))
491 gmx_ana_pos_append(NULL, out->u.g, pos, b, 0);
497 * \param[in] sel Selection element to query.
498 * \returns TRUE if the covered fraction can be estimated for \p sel with
499 * _gmx_selelem_estimate_coverfrac(), FALSE otherwise.
502 _gmx_selelem_can_estimate_cover(t_selelem *sel)
508 if (sel->type == SEL_BOOLEAN && sel->u.boolt == BOOL_OR)
517 if (child->type == SEL_EXPRESSION)
519 if (child->u.expr.method->name == sm_insolidangle.name)
521 if (bFound || bDynFound)
527 else if (child->u.expr.method
528 && (child->u.expr.method->flags & SMETH_DYNAMIC))
537 else if (!_gmx_selelem_can_estimate_cover(child))
547 * \param[in] sel Selection for which the fraction should be calculated.
548 * \returns Fraction of angles covered by the selection (between zero and one).
550 * The return value is undefined if _gmx_selelem_can_estimate_cover() returns
552 * Should be called after gmx_ana_evaluate_selections() has been called for the
556 _gmx_selelem_estimate_coverfrac(t_selelem *sel)
561 if (sel->type == SEL_EXPRESSION && sel->u.expr.method->name == sm_insolidangle.name)
563 t_methoddata_insolidangle *d = (t_methoddata_insolidangle *)sel->u.expr.mdata;
566 d->cfrac = estimate_covered_fraction(d);
570 if (sel->type == SEL_BOOLEAN && sel->u.boolt == BOOL_NOT)
572 cfrac = _gmx_selelem_estimate_coverfrac(sel->child);
580 /* Here, we assume that the selection is simple enough */
584 cfrac = _gmx_selelem_estimate_coverfrac(child);
595 * \param[in] x1 Unit vector 1.
596 * \param[in] x2 Unit vector 2.
597 * \returns Minus the dot product of \p x1 and \p x2.
599 * This function is used internally to calculate the distance between the
600 * unit vectors \p x1 and \p x2 to find out whether \p x2 is within the
601 * cone centered at \p x1. Currently, the cosine of the angle is used
602 * for efficiency, and the minus is there to make it behave like a normal
603 * distance (larger values mean longer distances).
606 sph_distc(rvec x1, rvec x2)
608 return -iprod(x1, x2);
612 * \param[in] p Partition to search.
613 * \param[in] value Value to search for.
614 * \returns The partition index in \p p that contains \p value.
616 * If \p value is outside the range of \p p, the first/last index is returned.
617 * Otherwise, the return value \c i satisfies \c p->p[i].left<=value and
618 * \c p->p[i+1].left>value
621 find_partition_bin(t_partition *p, real value)
623 int pmin, pmax, pbin;
625 /* Binary search the partition */
626 pmin = 0; pmax = p->n;
627 while (pmax > pmin + 1)
629 pbin = pmin + (pmax - pmin) / 2;
630 if (p->p[pbin].left <= value)
644 * \param[in] surf Surface data structure to search.
645 * \param[in] x Unit vector to find.
646 * \returns The bin index that contains \p x.
648 * The return value is an index to the \p surf->bin array.
651 find_surface_bin(t_methoddata_insolidangle *surf, rvec x)
657 phi = atan2(x[YY], x[XX]);
658 tbin = floor(theta / surf->tbinsize);
659 if (tbin >= surf->ntbins)
661 tbin = surf->ntbins - 1;
663 pbin = find_partition_bin(&surf->tbin[tbin], phi);
664 return surf->tbin[tbin].p[pbin].bin;
668 * \param[in,out] surf Surface data structure.
670 * Clears the reference points from the bins and (re)initializes the edges
671 * of the azimuthal bins.
674 clear_surface_points(t_methoddata_insolidangle *surf)
679 for (i = 0; i < surf->ntbins; ++i)
681 c = min(sin(surf->tbinsize*i), sin(surf->tbinsize*(i+1)))
682 * M_2PI / surf->targetbinsize + 1;
688 for (j = 0; j < c; ++j)
690 surf->tbin[i].p[j].left = -M_PI + j*M_2PI/c - 0.0001;
691 surf->tbin[i].p[j].bin = surf->nbins;
692 surf->bin[surf->nbins].n = 0;
695 surf->tbin[i].p[c].left = M_PI + 0.0001;
696 surf->tbin[i].p[c].bin = -1;
701 * \param[in,out] surf Surface data structure.
704 free_surface_points(t_methoddata_insolidangle *surf)
708 for (i = 0; i < surf->nbins; ++i)
712 sfree(surf->bin[i].x);
714 surf->bin[i].n_alloc = 0;
715 surf->bin[i].x = NULL;
720 * \param[in,out] surf Surface data structure.
721 * \param[in] tbin Bin number in the zenith angle direction.
722 * \param[in] pbin Bin number in the azimuthal angle direction.
723 * \param[in] x Point to store.
726 add_surface_point(t_methoddata_insolidangle *surf, int tbin, int pbin, rvec x)
730 bin = surf->tbin[tbin].p[pbin].bin;
731 /* Return if bin is already completely covered */
732 if (surf->bin[bin].n == -1)
734 /* Allocate more space if necessary */
735 if (surf->bin[bin].n == surf->bin[bin].n_alloc) {
736 surf->bin[bin].n_alloc += 10;
737 srenew(surf->bin[bin].x, surf->bin[bin].n_alloc);
739 /* Add the point to the bin */
740 copy_rvec(x, surf->bin[bin].x[surf->bin[bin].n]);
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.
750 mark_surface_covered(t_methoddata_insolidangle *surf, int tbin, int pbin)
754 bin = surf->tbin[tbin].p[pbin].bin;
755 surf->bin[bin].n = -1;
759 * \param[in,out] surf Surface data structure.
760 * \param[in] tbin Bin number in the zenith angle direction.
761 * \param[in] phi Azimuthal angle of \p x.
762 * \param[in] pdelta1 Width of the cone at the lower edge of \p tbin.
763 * \param[in] pdelta2 Width of the cone at the uppper edge of \p tbin.
764 * \param[in] pdeltamax Max. width of the cone inside \p tbin.
765 * \param[in] x Point to store (should have unit length).
768 update_surface_bin(t_methoddata_insolidangle *surf, int tbin,
769 real phi, real pdelta1, real pdelta2, real pdeltamax,
772 real pdelta, phi1, phi2;
773 int pbin1, pbin2, pbin;
775 /* Find the edges of the bins affected */
776 pdelta = max(max(pdelta1, pdelta2), pdeltamax);
787 pbin1 = find_partition_bin(&surf->tbin[tbin], phi1);
788 pbin2 = find_partition_bin(&surf->tbin[tbin], phi2);
789 /* Find the edges of completely covered region */
790 pdelta = min(pdelta1, pdelta2);
797 /* Loop over all affected bins */
801 /* Wrap bin around if end reached */
802 if (pbin == surf->tbin[tbin].n)
808 /* Check if bin is completely covered and update */
809 if (surf->tbin[tbin].p[pbin].left >= phi1
810 && surf->tbin[tbin].p[pbin+1].left <= phi2)
812 mark_surface_covered(surf, tbin, pbin);
816 add_surface_point(surf, tbin, pbin, x);
819 while (pbin++ != pbin2); /* Loop including pbin2 */
823 * \param[in,out] surf Surface data structure.
824 * \param[in] x Point to store (should have unit length).
826 * Finds all the bins covered by the cone centered at \p x and calls
827 * update_surface_bin() to update them.
830 store_surface_point(t_methoddata_insolidangle *surf, rvec x)
833 real pdeltamax, tmax;
834 real theta1, theta2, pdelta1, pdelta2;
838 phi = atan2(x[YY], x[XX]);
839 /* Find the maximum extent in the phi direction */
840 if (theta <= surf->angcut)
845 else if (theta >= M_PI - surf->angcut)
852 pdeltamax = asin(sin(surf->angcut) / sin(theta));
853 tmax = acos(cos(theta) / cos(surf->angcut));
855 /* Find the first affected bin */
856 tbin = max(floor((theta - surf->angcut) / surf->tbinsize), 0.0);
857 theta1 = tbin * surf->tbinsize;
858 if (theta1 < theta - surf->angcut)
866 /* Loop through all affected bins */
867 while (tbin < ceil((theta + surf->angcut) / surf->tbinsize)
868 && tbin < surf->ntbins)
870 /* Calculate the next boundaries */
871 theta2 = (tbin+1) * surf->tbinsize;
872 if (theta2 > theta + surf->angcut)
876 else if (tbin == surf->ntbins - 1)
882 pdelta2 = 2*asin(sqrt(
883 (sqr(sin(surf->angcut/2)) - sqr(sin((theta2-theta)/2))) /
884 (sin(theta) * sin(theta2))));
887 if (tmax >= theta1 && tmax <= theta2)
889 update_surface_bin(surf, tbin, phi, pdelta1, pdelta2, pdeltamax, x);
893 update_surface_bin(surf, tbin, phi, pdelta1, pdelta2, 0, x);
903 * \param[in,out] surf Surface data structure.
905 * Currently, this function does nothing.
908 optimize_surface_points(t_methoddata_insolidangle *surf)
910 /* TODO: Implement */
914 * \param[in] surf Surface data structure.
915 * \returns An estimate for the area covered by the reference points.
918 estimate_covered_fraction(t_methoddata_insolidangle *surf)
921 real cfrac, tfrac, pfrac;
924 for (t = 0; t < surf->ntbins; ++t)
926 tfrac = cos(t * surf->tbinsize) - cos((t+1) * surf->tbinsize);
927 for (p = 0; p < surf->tbin[t].n; ++p)
929 pfrac = surf->tbin[t].p[p+1].left - surf->tbin[t].p[p].left;
930 n = surf->bin[surf->tbin[t].p[p].bin].n;
931 if (n == -1) /* Bin completely covered */
933 cfrac += tfrac * pfrac;
935 else if (n > 0) /* Bin partially covered */
937 cfrac += tfrac * pfrac / 2; /* A rough estimate */
941 return cfrac / (4*M_PI);
945 * \param[in] surf Surface data structure to search.
946 * \param[in] x Unit vector to check.
947 * \returns TRUE if \p x is within the solid angle, FALSE otherwise.
950 is_surface_covered(t_methoddata_insolidangle *surf, rvec x)
954 bin = find_surface_bin(surf, x);
955 /* Check for completely covered bin */
956 if (surf->bin[bin].n == -1)
960 /* Check each point that partially covers the bin */
961 for (i = 0; i < surf->bin[bin].n; ++i)
963 if (sph_distc(x, surf->bin[bin].x[i]) < surf->distccut)